Initial commitHEADmaster

15 files changed, 2584 insertions, 0 deletions

diff --git a/src/Behavior.cpp b/src/Behavior.cpp
new file mode 100644
index 0000000..15b2ad9
--- /dev/null
+++ b/src/Behavior.cpp
@@ -0,0 +1,59 @@
+#include <HyperNeat/Behavior.hpp>
+#include <HyperNeat/NoveltyMetric.hpp>
+
+using namespace hyperneat;
+
+void
+Behavior::clear()
+{
+    _characterization.assign(_characterization.size(), 0.0);
+}
+
+void
+Behavior::reset(bool archive)
+{
+    if (_toBeArchived && archive) {
+        _noveltyMetric->_archive.emplace_back(_characterization);
+    }
+
+    _noveltyScore    = 0.0;
+    _toBeArchived    = false;
+    _criteriaReached = _noveltyMetric->_prms._criteriaReachedByDefault;
+    _characterization.assign(_characterization.size(), 0.0);
+}
+
+double&
+Behavior::at(size_t i)
+{
+    return _characterization[i];
+}
+
+double
+Behavior::getNoveltyScore() const
+{
+    return _noveltyScore;
+}
+
+bool
+Behavior::isToBeArchived() const
+{
+    return _toBeArchived;
+}
+
+Organism&
+Behavior::getOrganism()
+{
+    return *_organism;
+}
+
+const NoveltyMetric&
+Behavior::getNoveltyMetric() const
+{
+    return *_noveltyMetric;
+}
+
+const Vector<double>&
+Behavior::getCharacterization() const
+{
+    return _characterization;
+}
diff --git a/src/Cppn.cpp b/src/Cppn.cpp
new file mode 100644
index 0000000..fb1d19d
--- /dev/null
+++ b/src/Cppn.cpp
@@ -0,0 +1,292 @@
+#include <cmath>
+#include <HyperNeat/Cppn.hpp>
+#include <HyperNeat/Genome.hpp>
+#include <HyperNeat/QuadTree.hpp>
+#include <HyperNeat/Utils/Pi.hpp>
+#include <HyperNeat/NodeSearchPrms.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+void
+Cppn::create(const Genome& genome)
+{
+    using DepthSorter = Map<double, Map<size_t, const Genome::NodeGene*>>;
+    DepthSorter dSorter;
+
+    for (auto& i : genome._nodeGenes) {
+        dSorter[i.second._depth][i.first] = &i.second;
+    }
+
+    _inputs.resize(genome._inputs, 0.0);
+    _nodeLayers.resize(dSorter.size());
+
+    {
+        auto layerIter = _nodeLayers.begin();
+
+        for (auto& i : dSorter) {
+            layerIter->resize(i.second.size());
+            auto nodeIter = layerIter->begin();
+
+            for (auto& j : i.second) {
+                size_t enabledLinks = 0;
+
+                for (auto& k : j.second->_linkGenes) {
+                    if (k.second._isEnabled) {
+                        ++enabledLinks;
+                    }
+                }
+
+                nodeIter->_links.resize(enabledLinks);
+                ++nodeIter;
+            }
+
+            ++layerIter;
+        }
+    }
+
+    size_t layerIdx = 0;
+
+    for (auto& i : dSorter) {
+        size_t nodeIdx = 0;
+
+        for (auto& j : i.second) {
+            Node& crntNode = _nodeLayers[layerIdx][nodeIdx];
+            size_t linkIdx = 0;
+
+            for (auto& k : j.second->_linkGenes) {
+                if (k.second._isEnabled) {
+                    Node::Link& crntLink = crntNode._links[linkIdx];
+
+                    if (k.first < _inputs.size()) {
+                        crntLink._input = &_inputs[k.first];
+                    } else {
+                        double depth = genome._nodeGenes.at(k.first)._depth;
+                        size_t layer = distance(dSorter.begin(), dSorter.find(depth));
+                        size_t node  = distance(dSorter[depth].begin(), dSorter[depth].find(k.first));
+
+                        crntLink._input = &_nodeLayers[layer][node]._output;
+                    }
+
+                    crntLink._weight = k.second._weight;
+                    ++linkIdx;
+                }
+            }
+
+            crntNode._nodeType = j.second->_nodeType;
+            ++nodeIdx;
+
+            if (i.first == 1.0) {
+                _outputs.push_back(&crntNode._output);
+            }
+        }
+
+        ++layerIdx;
+    }
+}
+
+void
+Cppn::clear()
+{
+    _inputs.clear();
+    _outputs.clear();
+    _nodeLayers.clear();
+}
+
+size_t
+Cppn::getInputsCount() const
+{
+    return _inputs.size();
+}
+
+size_t
+Cppn::getOutputsCount() const
+{
+    return _outputs.size();
+}
+
+size_t
+Cppn::getNodesCount() const
+{
+    size_t nodesCount = 0;
+
+    for (auto& i : _nodeLayers) {
+        nodesCount += i.size();
+    }
+
+    return nodesCount;
+}
+
+double&
+Cppn::inputAt(size_t i)
+{
+    return _inputs[i];
+}
+
+double
+Cppn::outputAt(size_t i) const
+{
+    return *_outputs[i];
+}
+
+void
+Cppn::cycle()
+{
+    for (auto& i : _nodeLayers) {
+        for (auto& j : i) {
+            j.appendInput();
+        }
+
+        for (auto& j : i) {
+            j.flushOutput();
+        }
+    }
+}
+
+void
+Cppn::findNodesIn2DSection(ValueMap& valueMap, const NodeSearchPrms& nsPrms, const Point& source)
+{
+    QuadTree qTree(1.0, 0.0, 0.0);
+
+    qTree.subdivide([&](QuadTree* qt) {
+        auto levelToSegment = [](size_t level) {
+            double dlevel = static_cast<double>(level);
+            return 1.0 / pow(2.0, dlevel);
+        };
+
+        double minQuadTreeSegment = levelToSegment(nsPrms._maxQuadTreeLevel);
+        double maxQuadTreeSegment = levelToSegment(nsPrms._minQuadTreeLevel);
+        double testGridResolution = levelToSegment(nsPrms._testGridLevel);
+
+        if (qt->getSegment() <= minQuadTreeSegment) {
+            return false;
+        }
+
+        if (qt->getSegment() > maxQuadTreeSegment) {
+            return true;
+        }
+
+        double gridStep = qt->getSegment() * 2.0 * testGridResolution;
+        double gridHalf = gridStep / 2.0;
+        double initX    = qt->getX() - qt->getSegment() + gridHalf;
+        double initY    = qt->getY() - qt->getSegment() + gridHalf;
+        double endX     = qt->getX() + qt->getSegment();
+        double endY     = qt->getY() + qt->getSegment();
+        double cellsNo  = (1.0 / testGridResolution) * (1.0 / testGridResolution);
+
+        Vector<double> w;
+        w.reserve(static_cast<size_t>(cellsNo));
+
+        for (double y = initY; y < endY; y += gridStep) {
+            for (double x = initX; x < endX; x += gridStep) {
+                _inputs[nsPrms._x] = x;
+                _inputs[nsPrms._y] = y;
+
+                if (nsPrms._useDistance) {
+                    _inputs[nsPrms._d] = Point(x, y).distance(source);
+                }
+
+                cycle();
+                w.emplace_back(*_outputs[nsPrms._o]);
+            }
+        }
+
+        double wMean    = 0.0;
+        double variance = 0.0;
+
+        for (auto& i : w) {
+            wMean += i;
+        }
+
+        wMean /= cellsNo;
+
+        for (auto& i : w) {
+            variance += (wMean - i) * (wMean - i);
+        }
+
+        variance /= cellsNo;
+
+        if (variance > nsPrms._varianceThreshold) {
+            return true;
+        }
+
+        return false;
+    });
+
+    qTree.traverse([&](const QuadTree* qt) {
+        auto measureDistAndIO = [&]() {
+            if (nsPrms._useDistance) {
+                _inputs[nsPrms._d] = Point(_inputs[nsPrms._x], _inputs[nsPrms._y]).distance(source);
+            }
+
+            cycle();
+            return *_outputs[nsPrms._o];
+        };
+
+        double gridSize  = qt->getSegment() * 2.0;
+
+        _inputs[nsPrms._y] = qt->getY();
+        _inputs[nsPrms._x] = qt->getX() - gridSize;
+        double valLeft     = measureDistAndIO();
+
+        _inputs[nsPrms._x] = qt->getX() + gridSize;
+        double valRight    = measureDistAndIO();
+
+        _inputs[nsPrms._x] = qt->getX();
+        _inputs[nsPrms._y] = qt->getY() - gridSize;
+        double valBottom   = measureDistAndIO();
+
+        _inputs[nsPrms._y] = qt->getY() + gridSize;
+        double valTop      = measureDistAndIO();
+
+        _inputs[nsPrms._y] = qt->getY();
+        double valCenter   = measureDistAndIO();
+
+        valLeft   = fabs(valCenter - valLeft);
+        valRight  = fabs(valCenter - valRight);
+        valBottom = fabs(valCenter - valBottom);
+        valTop    = fabs(valCenter - valTop);
+
+        double bandValue = max(min(valLeft, valRight), min(valBottom, valTop));
+
+        if (bandValue > nsPrms._bandPruningThreshold) {
+            valueMap.emplace_back(qt->getX(), qt->getY(), valCenter, qt->getSegment());
+        }
+    });
+}
+
+void
+Cppn::Node::appendInput()
+{
+    for (auto& i : _links) {
+        _storedInput += *i._input * i._weight;
+    }
+}
+
+void
+Cppn::Node::flushOutput()
+{
+    switch (_nodeType) {
+        case NodeType::SIGMOID:
+            _output = 1.0 / (1.0 + exp(-_storedInput * 2.0));
+            _output = _output * 2.0 - 1.0;
+            break;
+
+        case NodeType::GAUSSIAN:
+            _output = exp(-_storedInput * _storedInput);
+            _output = _output * 2.0 - 1.0;
+            break;
+
+        case NodeType::SINE:
+            _output = sin((_storedInput / 2.0) * PI);
+            break;
+
+        case NodeType::ABSOLUTE:
+            _output = max(min(fabs(_storedInput), 0.0), 1.0);
+            break;
+
+        default:;
+    }
+
+    _storedInput = 0.0;
+}
diff --git a/src/Genome.cpp b/src/Genome.cpp
new file mode 100644
index 0000000..ce5fc11
--- /dev/null
+++ b/src/Genome.cpp
@@ -0,0 +1,15 @@
+#include <HyperNeat/Genome.hpp>
+
+using namespace hyperneat;
+
+Genome::Genome(size_t inputs)
+    : _inputs(inputs)
+{}
+
+Genome::NodeGene::NodeGene(double depth, NodeType nodeType)
+    : _depth(depth), _nodeType(nodeType)
+{}
+
+Genome::NodeGene::LinkGene::LinkGene(double weight, bool isEnabled)
+    : _weight(weight), _isEnabled(isEnabled)
+{}
diff --git a/src/Innovation.cpp b/src/Innovation.cpp
new file mode 100644
index 0000000..0332616
--- /dev/null
+++ b/src/Innovation.cpp
@@ -0,0 +1,16 @@
+#include <Hyperneat/Innovation.hpp>
+
+using namespace hyperneat;
+
+Innovation::Innovation(size_t number, size_t source, size_t target, double depth, NodeType nodeType)
+    : _number(number), _source(source), _target(target), _depth(depth), _nodeType(nodeType)
+{}
+
+bool
+Innovation::operator==(const Innovation& other) const
+{
+    return (_source   == other._source  ) &&
+           (_target   == other._target  ) &&
+           (_depth    == other._depth   ) &&
+           (_nodeType == other._nodeType);
+}
diff --git a/src/NeuralNet.cpp b/src/NeuralNet.cpp
new file mode 100644
index 0000000..abed262
--- /dev/null
+++ b/src/NeuralNet.cpp
@@ -0,0 +1,294 @@
+#include <cmath>
+#include <algorithm>
+#include <HyperNeat/Cppn.hpp>
+#include <HyperNeat/NeuralNet.hpp>
+#include <HyperNeat/Utils/Map.hpp>
+#include <HyperNeat/Utils/Set.hpp>
+#include <HyperNeat/NeuralNetPrms.hpp>
+#include <HyperNeat/Utils/Function.hpp>
+#include <HyperNeat/NodeSearchPrms.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+void
+NeuralNet::create(Cppn& cppn, const NeuralNetPrms& nnPrms)
+{
+    class TempNeuronSynapse {
+    public:
+        TempNeuronSynapse() = default;
+        TempNeuronSynapse(double weight, Point neuronSource)
+            : _weight(weight), _neuronSource(neuronSource)
+        {}
+
+        double _weight = 0.0;
+        Point  _neuronSource;
+    };
+
+    class TempNeuron {
+    public:
+        bool                      _isIncluded = false;
+        double                    _bias       = 0.0;
+        Vector<TempNeuronSynapse> _neuronSynapses;
+    };
+
+    _inputs.reserve(nnPrms._inputMap.size());
+    _outputs.reserve(nnPrms._outputMap.size());
+
+    NodeSearchPrms nsPrms(0, 2, 3, 4);
+    nsPrms.importFrom(nnPrms);
+    cppn.inputAt(5) = 1.0;
+
+    Map<Point, TempNeuron> tempNeurons;
+
+    auto findConnections = [&](const Point& source, size_t x1, size_t y1, size_t x2, size_t y2, size_t d,
+            bool checkExist, Function<void(double, const Point&)> storeConn) {
+        cppn.inputAt(x1) = source._x;
+        cppn.inputAt(y1) = source._y;
+        cppn.inputAt(x2) = 0.0;
+        cppn.inputAt(y2) = 0.0;
+        cppn.inputAt(d)  = 0.0;
+
+        ValueMap newConnections;
+        cppn.findNodesIn2DSection(newConnections, nsPrms, source);
+
+        for (auto& i : newConnections) {
+            if (checkExist && !tempNeurons.count(i)) {
+                continue;
+            }
+
+            if (fabs(i._value) > 0.2) {
+                storeConn((i._value + (i._value > 0 ? -0.2 : 0.2)) * 3.75, i);
+            }
+        }
+    };
+
+    {
+        Set<Point>  neuronSet1;
+        Set<Point>  neuronSet2;
+        Set<Point>* previousNeurons = &neuronSet1;
+        Set<Point>* nextNeurons     = &neuronSet2;
+
+        for (auto& i : nnPrms._inputMap) {
+            tempNeurons[i];
+            previousNeurons->insert(i);
+        }
+
+        for (size_t i = 0; i < nnPrms._iterations && !previousNeurons->empty(); ++i) {
+            Map<Point, TempNeuron> newNeurons;
+
+            for (auto& j : *previousNeurons) {
+                findConnections(j, 0, 1, 2, 3, 4, false, [&](double weight, const Point& target) {
+                    if (tempNeurons.count(target)) {
+                        auto& synapses = tempNeurons[target]._neuronSynapses;
+                        synapses.emplace_back(weight, j);
+                    } else {
+                        auto& synapses = newNeurons[target]._neuronSynapses;
+                        synapses.emplace_back(weight, j);
+                        nextNeurons->insert(target);
+                    }
+                });
+            }
+
+            previousNeurons->clear();
+            swap(nextNeurons, previousNeurons);
+            tempNeurons.insert(newNeurons.begin(), newNeurons.end());
+        }
+    }
+
+    nsPrms._x = 0;
+    nsPrms._y = 1;
+
+    {
+        Vector<TempNeuron*>  inclusionSet1;
+        Vector<TempNeuron*>  inclusionSet2;
+        Vector<TempNeuron*>* crntInclusions = &inclusionSet1;
+        Vector<TempNeuron*>* nextInclusions = &inclusionSet2;
+
+        for (auto& i : nnPrms._outputMap) {
+            tempNeurons[i]._isIncluded = true;
+            nextInclusions->push_back(&tempNeurons[i]);
+
+            findConnections(i, 2, 3, 0, 1, 4, true, [&](double weight, const Point& target) {
+                auto& synapses = tempNeurons[i]._neuronSynapses;
+                synapses.emplace_back(weight, target);
+            });
+        }
+
+        while (!nextInclusions->empty()) {
+            crntInclusions->clear();
+            swap(crntInclusions, nextInclusions);
+
+            for (auto& i : *crntInclusions) {
+                for (auto& j : i->_neuronSynapses) {
+                    auto& sourceNeuron = tempNeurons.at(j._neuronSource);
+
+                    if (!sourceNeuron._isIncluded) {
+                        nextInclusions->push_back(&sourceNeuron);
+                        sourceNeuron._isIncluded = true;
+                    }
+                }
+            }
+        }
+    }
+
+    for (auto& i : nnPrms._inputMap) {
+        tempNeurons[i]._isIncluded = true;
+    }
+
+    cppn.inputAt(2) = 0.0;
+    cppn.inputAt(3) = 0.0;
+    cppn.inputAt(4) = 0.0;
+
+    for (auto i = tempNeurons.begin(), end = tempNeurons.end(); i != end;) {
+        if (i->second._isIncluded) {
+            cppn.inputAt(0) = i->first._x;
+            cppn.inputAt(1) = i->first._y;
+            cppn.inputAt(4) = i->first.distance(Point());
+            cppn.cycle();
+            i->second._bias = cppn.outputAt(1) * 3.0;
+            ++i;
+        } else {
+            i = tempNeurons.erase(i);
+        }
+    }
+
+    _neurons.resize(tempNeurons.size());
+
+    {
+        auto nIter = _neurons.begin();
+
+        for (auto& i : tempNeurons) {
+            nIter->_bias         = i.second._bias;
+            nIter->_position     = i.first;
+            auto& crntNrnSyns    = nIter->_synapses;
+            auto& neuronSynapses = i.second._neuronSynapses;
+            crntNrnSyns.reserve(neuronSynapses.size());
+
+            for (auto& j : neuronSynapses) {
+                auto src    = tempNeurons.find(j._neuronSource);
+                size_t sIdx = distance(tempNeurons.begin(), src);
+                crntNrnSyns.emplace_back(&_neurons[sIdx], j._weight);
+            }
+
+            ++nIter;
+        }
+    }
+
+    auto relateIO = [&](Vector<double*>& ptrVec, const Vector<Point>& map, Neuron::Type type) {
+        for (auto& i : map) {
+            auto neuron          = tempNeurons.find(i);
+            size_t nIdx          = distance(tempNeurons.begin(), neuron);
+            _neurons[nIdx]._type = type;
+            ptrVec.push_back(&_neurons[nIdx]._output);
+        }
+    };
+
+    relateIO(_inputs, nnPrms._inputMap, Neuron::Type::INPUT);
+    relateIO(_outputs, nnPrms._outputMap, Neuron::Type::OUTPUT);
+}
+
+void
+NeuralNet::clear()
+{
+    _inputs.clear();
+    _outputs.clear();
+    _neurons.clear();
+}
+
+size_t
+NeuralNet::getInputsCount() const
+{
+    return _inputs.size();
+}
+
+size_t
+NeuralNet::getOutputsCount() const
+{
+    return _outputs.size();
+}
+
+size_t
+NeuralNet::getNeuronsCount() const
+{
+    return _neurons.size();
+}
+
+double
+NeuralNet::getAverageActivation() const
+{
+    double totalActivation = 0.0;
+
+    for (auto& i : _neurons) {
+        totalActivation += i._output;
+    }
+
+    return totalActivation / static_cast<double>(_neurons.size());
+}
+
+double&
+NeuralNet::inputAt(size_t i)
+{
+    return *_inputs[i];
+}
+
+double
+NeuralNet::outputAt(size_t i) const
+{
+    return *_outputs[i];
+}
+
+const Vector<double*>&
+NeuralNet::getInputs() const
+{
+    return _inputs;
+}
+
+const Vector<double*>&
+NeuralNet::getOutputs() const
+{
+    return _outputs;
+}
+
+const Vector<NeuralNet::Neuron>&
+NeuralNet::getNeurons() const
+{
+    return _neurons;
+}
+
+void
+NeuralNet::cycle()
+{
+    for (auto& i : _neurons) {
+        i.appendInput();
+    }
+
+    for (auto& i : _neurons) {
+        i.flushOutput();
+    }
+}
+
+NeuralNet::Neuron::Neuron(const Point& position, Type type, double bias)
+    : _position(position), _type(type), _bias(bias)
+{}
+
+void
+NeuralNet::Neuron::appendInput()
+{
+    for (auto& i : _synapses) {
+        _storedInput += *i._input * i._weight;
+    }
+
+    _storedInput += _bias;
+}
+
+void
+NeuralNet::Neuron::flushOutput()
+{
+    _output      = 1.0 / (1.0 + exp(-_storedInput * 4.9));
+    _storedInput = 0.0;
+}
+
+NeuralNet::Neuron::Synapse::Synapse(Neuron* inputNeuron, double weight)
+    : _input(&inputNeuron->_output), _neuron(inputNeuron), _weight(weight)
+{}
diff --git a/src/NodeSearchPrms.cpp b/src/NodeSearchPrms.cpp
new file mode 100644
index 0000000..db2421a
--- /dev/null
+++ b/src/NodeSearchPrms.cpp
@@ -0,0 +1,23 @@
+#include <HyperNeat/NeuralNetPrms.hpp>
+#include <HyperNeat/NodeSearchPrms.hpp>
+
+using namespace hyperneat;
+
+NodeSearchPrms::NodeSearchPrms(size_t o, size_t x, size_t y)
+    : _o(o), _x(x), _y(y), _useDistance(false)
+{}
+
+NodeSearchPrms::NodeSearchPrms(size_t o, size_t x, size_t y, size_t d)
+    : _o(o), _x(x), _y(y), _d(d)
+{}
+
+void
+NodeSearchPrms::importFrom(const NeuralNetPrms& nnPrms)
+{
+    _testGridLevel        = nnPrms._testGridLevel;
+    _maxQuadTreeLevel     = nnPrms._maxQuadTreeLevel;
+    _minQuadTreeLevel     = nnPrms._minQuadTreeLevel;
+    _bandPruningThreshold = nnPrms._bandPruningThreshold;
+    _varianceThreshold    = nnPrms._varianceThreshold;
+    _divisionThreshold    = nnPrms._divisionThreshold;
+}
diff --git a/src/NoveltyMetric.cpp b/src/NoveltyMetric.cpp
new file mode 100644
index 0000000..dd40118
--- /dev/null
+++ b/src/NoveltyMetric.cpp
@@ -0,0 +1,101 @@
+#include <algorithm>
+#include <Hyperneat/Population.hpp>
+#include <Hyperneat/NoveltyMetric.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+const NoveltyMetricPrms&
+NoveltyMetric::getPrms() const
+{
+    return _prms;
+}
+
+const Vector<Behavior>&
+NoveltyMetric::getBehaviors() const
+{
+    return _behaviors;
+}
+
+const Vector2D<double>&
+NoveltyMetric::getArchive() const
+{
+    return _archive;
+}
+
+Behavior&
+NoveltyMetric::getBehaviorOf(size_t i)
+{
+    return _behaviors[i];
+}
+
+void
+NoveltyMetric::initialize(const NoveltyMetricPrms& prms, Population* population)
+{
+    _prms = prms;
+    _behaviors.resize(population->getAllOrganisms().size());
+
+    for (size_t i = 0; i < _behaviors.size(); ++i) {
+        population->getOrganism(i)._behavior = &_behaviors[i];
+        _behaviors[i]._organism              = &population->getOrganism(i);
+        _behaviors[i]._noveltyMetric         = this;
+        _behaviors[i]._criteriaReached       = prms._criteriaReachedByDefault;
+        _behaviors[i]._characterization.resize(prms._characterizationSize, 0.0);
+    }
+}
+
+void
+NoveltyMetric::setScores()
+{
+    size_t readyOrganisms = 0;
+
+    for (auto& i : _behaviors) {
+        if ((i._isReady = (i._organism->isOld() && i._criteriaReached))) {
+            ++readyOrganisms;
+        }
+    }
+
+    for (auto& i : _behaviors) {
+        if (!i._isReady) {
+            i._toBeArchived       = false;
+            i._organism->_fitness = i._noveltyScore = 0.0;
+            continue;
+        }
+
+        Vector<double> distances;
+        distances.reserve(readyOrganisms + _archive.size() - 1);
+
+        for (auto& j : _behaviors) {
+            if (j._isReady && (&i != &j)) {
+                distances.emplace_back(getDistance(i._characterization, j._characterization));
+            }
+        }
+
+        for (auto& j : _archive) {
+            distances.emplace_back(getDistance(i._characterization, j));
+        }
+
+        auto dEnd = distances.end();
+        auto dBeg = distances.begin();
+        auto dRef = dBeg + _prms._referenceOrganisms;
+
+        partial_sort(dBeg, dRef, dEnd);
+
+        double totalDistance  = accumulate(dBeg, dRef, 0.0);
+        i._organism->_fitness = i._noveltyScore = (totalDistance / static_cast<double>(_prms._referenceOrganisms));
+        i._toBeArchived       = (i._noveltyScore > _prms._noveltyThreshold);
+    }
+}
+
+double
+NoveltyMetric::getDistance(const Vector<double>& v1, const Vector<double>& v2) const
+{
+    double result = 0.0;
+
+    for (size_t i = 0; i < v1.size(); ++i) {
+        double diff = v1[i] - v2[i];
+        result     += diff * diff;
+    }
+
+    return sqrt(result);
+}
diff --git a/src/Organism.cpp b/src/Organism.cpp
new file mode 100644
index 0000000..0c2c1bc
--- /dev/null
+++ b/src/Organism.cpp
@@ -0,0 +1,180 @@
+#include <HyperNeat/Cppn.hpp>
+#include <HyperNeat/Behavior.hpp>
+#include <HyperNeat/NeuralNet.hpp>
+#include <HyperNeat/Population.hpp>
+#include <HyperNeat/Utils/Thread.hpp>
+#include <HyperNeat/NeuralNetPrms.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+Organism::Organism(const Organism& other)
+{
+    *this = other;
+}
+
+Organism&
+Organism::operator=(const Organism& other)
+{
+    _fitness    = other._fitness;
+    _index      = other._index;
+    _isLocked   = other._isLocked;
+    _isFrozen   = other._isFrozen;
+    _specie     = other._specie;
+    _lifetime   = other._lifetime;
+    _genome     = other._genome;
+    _population = other._population;
+
+    return *this;
+}
+
+size_t
+Organism::getIndex() const
+{
+    return _index;
+}
+
+void
+Organism::lock()
+{
+    if (!_isLocked) {
+        _isLocked = true;
+        ++_population->_lockedOrganisms;
+    }
+}
+
+void
+Organism::unlock()
+{
+    if (_isLocked) {
+        _isLocked = false;
+        --_population->_lockedOrganisms;
+    }
+}
+
+bool
+Organism::isLocked() const
+{
+    return _isLocked;
+}
+
+void
+Organism::freeze()
+{
+    if (!_isFrozen) {
+        _isFrozen = true;
+        ++_population->_frozenOrganisms;
+    }
+}
+
+void
+Organism::unfreeze()
+{
+    if (_isFrozen) {
+        _isFrozen = false;
+        --_population->_frozenOrganisms;
+    }
+}
+
+bool
+Organism::isFrozen() const
+{
+    return _isFrozen;
+}
+
+bool
+Organism::isBeingGenerated() const
+{
+    return _isBeingGenerated;
+}
+
+size_t
+Organism::getSpecie() const
+{
+    return _specie;
+}
+
+bool
+Organism::isOld() const {
+    return _lifetime >= _population->_prms._minimumLifetime;
+}
+
+size_t
+Organism::getLifetime() const
+{
+    return _lifetime;
+}
+
+Behavior&
+Organism::getBehavior()
+{
+    return *_behavior;
+}
+
+const Genome&
+Organism::getGenome() const
+{
+    return _genome;
+}
+
+bool
+Organism::isChampion() const
+{
+    return &_population->getChampion() == this;
+}
+
+Population&
+Organism::getPopulation() const
+{
+    return *_population;
+}
+
+void
+Organism::createNeuralNet()
+{
+    if (_isBeingGenerated) {
+        return;
+    }
+
+    _isBeingGenerated = true;
+    ++_population->_organismsBeingGenerated;
+
+    Thread generator([&]() {
+        Cppn cppn;
+        cppn.create(_genome);
+
+        _neuralNet->clear();
+        _neuralNet->create(cppn, _population->getNeuralNetPrms());
+
+        _isBeingGenerated = false;
+        --_population->_organismsBeingGenerated;
+    });
+
+    generator.detach();
+}
+
+Organism::Organism(Population* population)
+    : _population(population)
+{}
+
+Organism::Organism(size_t inputs, Population* population)
+    : _genome(inputs), _population(population)
+{}
+
+void
+Organism::reset(bool archive)
+{
+    unlock();
+    unfreeze();
+
+    if (isOld()) {
+        --_population->_oldOrganisms;
+    }
+
+    _lifetime = 0;
+    _fitness  = 0.0;
+
+    if (_behavior) {
+        _behavior->reset(archive);
+    }
+}
diff --git a/src/Population.cpp b/src/Population.cpp
new file mode 100644
index 0000000..f81887e
--- /dev/null
+++ b/src/Population.cpp
@@ -0,0 +1,939 @@
+#include <chrono>
+#include <algorithm>
+#include <HyperNeat/Cppn.hpp>
+#include <HyperNeat/NeuralNet.hpp>
+#include <HyperNeat/Population.hpp>
+#include <HyperNeat/Utils/Thread.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+void
+Population::create(const PopulationPrms& popPrms)
+{
+    if (_prms._seed != 0) {
+        _randGen.seed(_prms._seed);
+    } else {
+        _randGen.seed(getRandSeed());
+    }
+
+    _prms              = popPrms;
+    _weightDeviator    = BellDist(0.0, _prms._weightDeviation);
+    _weightSelector    = RealDist(-_prms._weightRange, _prms._weightRange);
+    _distanceThreshold = _prms._initialDistanceThreshold;
+    _minOldOrganisms   = _prms._popSize * _prms._eligibilityRatio;
+    _basicInnovs       = _prms._cppnOutputs * NODE_TYPES_COUNT + _prms._cppnInputs;
+
+    _allOrganisms.resize(_prms._popSize, Organism(_prms._cppnInputs, this));
+
+    size_t orgIndex = 0;
+
+    for (auto& i : _allOrganisms) {
+        i._index        = orgIndex++;
+        auto& nodeGenes = i._genome._nodeGenes;
+
+        for (size_t j = 0; j < _prms._cppnOutputs; ++j) {
+            size_t nodeType    = _nodeTypeSelector(_randGen);
+            size_t geneIdx     = _prms._cppnInputs + (j * NODE_TYPES_COUNT) + nodeType;
+            nodeGenes[geneIdx] = {1.0, static_cast<NodeType>(nodeType)};
+
+            for (size_t k = 0; k < _prms._cppnInputs; ++k) {
+                nodeGenes[geneIdx]._linkGenes[k] = {getRandWeight()};
+            }
+        }
+    }
+
+    organizeSpecies();
+}
+
+void
+Population::create(const PopulationPrms& popPrms, const NeuralNetPrms& nnPrms)
+{
+    create(popPrms);
+    _nnPrms = Pointer<NeuralNetPrms>(new NeuralNetPrms(nnPrms));
+    generateAllNeuralNets();
+}
+
+void
+Population::create(const PopulationPrms& popPrms, const NeuralNetPrms& nnPrms, const NoveltyMetricPrms& nmPrms)
+{
+    create(popPrms, nnPrms);
+    setNoveltyMetric(nmPrms);
+}
+
+void
+Population::shutdown(bool resetOrganisms, bool archiveOrganisms)
+{
+    while (isAnyOrganismBeingGenerated());
+
+    if (resetOrganisms) {
+        for (auto& i : _allOrganisms) {
+            i.reset(archiveOrganisms);
+        }
+    }
+}
+
+Population::~Population()
+{
+    while (isAnyOrganismBeingGenerated());
+}
+
+void
+Population::setMinimumLifetime(size_t lifetime)
+{
+    _oldOrganisms          = 0;
+    _prms._minimumLifetime = lifetime;
+
+    for (auto& i : _allOrganisms) {
+        if (i._lifetime >= lifetime) {
+            ++_oldOrganisms;
+        }
+    }
+}
+
+const PopulationPrms&
+Population::getPopulationPrms() const
+{
+    return _prms;
+}
+
+const NeuralNetPrms&
+Population::getNeuralNetPrms() const
+{
+    return *_nnPrms;
+}
+
+bool
+Population::hasNeuralNets() const
+{
+    return (bool)_nnPrms;
+}
+
+const Vector<Organism>&
+Population::getAllOrganisms() const
+{
+    return _allOrganisms;
+}
+
+const Vector2D<Organism*>&
+Population::getSpecies() const
+{
+    return _species;
+}
+
+Organism&
+Population::getOrganism(size_t i)
+{
+    return _allOrganisms[i];
+}
+
+const Vector<Organism*>&
+Population::getSpecie(size_t i) const
+{
+    return _species[i];
+}
+
+Organism&
+Population::getChampion()
+{
+    auto comFitness = [](Organism& a, Organism& b) {
+        return a._fitness < b._fitness;
+    };
+
+    auto comNovelty = [](Organism& a, Organism& b) {
+        return a.getBehavior().getNoveltyScore() < b.getBehavior().getNoveltyScore();
+    };
+
+    return *max_element(_allOrganisms.begin(), _allOrganisms.end(), isNoveltyMetricSet() ? comNovelty : comFitness);
+}
+
+void
+Population::lock()
+{
+    _populationLock = true;
+}
+
+void
+Population::unlock()
+{
+    _populationLock = false;
+}
+
+bool
+Population::isLocked() const
+{
+    return _populationLock;
+}
+
+void
+Population::lockOrganism(size_t i)
+{
+    _allOrganisms[i].lock();
+}
+
+void
+Population::unlockOrganism(size_t i)
+{
+    _allOrganisms[i].unlock();
+}
+
+bool
+Population::isOrganismLocked(size_t i) const
+{
+    return _allOrganisms[i]._isLocked;
+}
+
+bool
+Population::isAnyOrganismLocked() const
+{
+    return _lockedOrganisms != 0 ? true : false;
+}
+
+size_t
+Population::getLockedOrganisms() const
+{
+    return _lockedOrganisms;
+}
+
+void
+Population::freezeOrganism(size_t i)
+{
+    _allOrganisms[i].freeze();
+}
+
+void
+Population::unfreezeOrganism(size_t i)
+{
+    _allOrganisms[i].unfreeze();
+}
+
+bool
+Population::isOrganismFrozen(size_t i) const
+{
+    return _allOrganisms[i]._isFrozen;
+}
+
+bool
+Population::isAnyOrganismFrozen() const
+{
+    return _frozenOrganisms != 0 ? true : false;
+}
+
+size_t
+Population::getFrozenOrganisms() const
+{
+    return _frozenOrganisms;
+}
+
+bool
+Population::isOrganismBeingGenerated(size_t i) const
+{
+    return _allOrganisms[i].isBeingGenerated();
+}
+
+bool
+Population::isAnyOrganismBeingGenerated() const
+{
+    return _organismsBeingGenerated != 0 ? true : false;
+}
+
+size_t
+Population::getOrganismsBeingGenerated() const
+{
+    return _organismsBeingGenerated;
+}
+
+size_t
+Population::getReadyOrganisms() const
+{
+    size_t result = 0;
+
+    for (auto& i : _allOrganisms) {
+        if (i.isOld() && !i.isLocked()) {
+            ++result;
+        }
+    }
+
+    return result;
+}
+
+const Vector<Innovation>&
+Population::getInnovations() const
+{
+    return _innovations;
+}
+
+size_t
+Population::getInnovationsCount() const
+{
+    return _innovations.size();
+}
+
+size_t
+Population::getBasicInnovationsCount() const
+{
+    return _basicInnovs;
+}
+
+Organism*
+Population::getLastReplacement()
+{
+    return _lastReplacement;
+}
+
+Organism*
+Population::getLastMother()
+{
+    return _lastMother;
+}
+
+Organism*
+Population::getLastFather()
+{
+    return _lastFather;
+}
+
+size_t
+Population::getReplacements() const
+{
+    return _replacements;
+}
+
+bool
+Population::recentReplacement() const
+{
+    return _recentReplacement;
+}
+
+double
+Population::getDistanceThreshold() const
+{
+    return _distanceThreshold;
+}
+
+size_t
+Population::getOldOrganisms() const
+{
+    return _oldOrganisms;
+}
+
+size_t
+Population::getMinimumOldOrganisms() const
+{
+    return _minOldOrganisms;
+}
+
+double
+Population::getAverageFitness() const
+{
+    double total = 0.0;
+
+    for (auto& i : _allOrganisms) {
+        total += i._fitness;
+    }
+
+    return total / static_cast<double>(_allOrganisms.size());
+}
+
+double
+Population::getAverageOldFitness() const
+{
+    double total = 0.0;
+
+    for (auto& i : _allOrganisms) {
+        if (i.isOld()) {
+            total += i._fitness;
+        }
+    }
+
+    return total / static_cast<double>(getOldOrganisms());
+}
+
+void
+Population::setNoveltyMetric(const NoveltyMetricPrms& prms) {
+    _noveltyMetric = Pointer<NoveltyMetric>(new NoveltyMetric);
+    _noveltyMetric->initialize(prms, this);
+}
+
+void
+Population::clearNoveltyMetric()
+{
+    _noveltyMetric.reset();
+
+    for (auto& i : _allOrganisms) {
+        i._behavior = nullptr;
+    }
+}
+
+bool
+Population::isNoveltyMetricSet() const
+{
+    return _noveltyMetric.get();
+}
+
+const NoveltyMetric&
+Population::getNoveltyMetric() const
+{
+    return *_noveltyMetric;
+}
+
+size_t
+Population::getUpdates() const
+{
+    return _updates;
+}
+
+double&
+Population::fitnessOf(size_t i)
+{
+    return _allOrganisms[i]._fitness;
+}
+
+bool
+Population::update(Function<void(void)> beforeReplacement, Function<void(void)> afterReplacement)
+{
+    ++_updates;
+
+    for (auto& i : _allOrganisms) {
+        if (!i._isFrozen && !i.isBeingGenerated()) {
+            ++i._lifetime;
+
+            if (i._lifetime == _prms._minimumLifetime) {
+                ++_oldOrganisms;
+            }
+        }
+    }
+
+    if (_populationLock) {
+        return false;
+    }
+
+    if ((getReadyOrganisms() >= _minOldOrganisms) && !isAnyOrganismBeingGenerated()) {
+        beforeReplacement();
+
+        if (isNoveltyMetricSet()) {
+            _noveltyMetric->setScores();
+        }
+
+        _recentReplacement = true;
+        replaceOrganism();
+        afterReplacement();
+    } else {
+        _recentReplacement = false;
+    }
+
+    return _recentReplacement;
+}
+
+void
+Population::generateAllNeuralNets()
+{
+    auto generateChunk = [&](size_t i, size_t end) {
+        for (; i < end; ++i) {
+            Cppn cppn;
+            cppn.create(_allOrganisms[i].getGenome());
+
+            _allOrganisms[i]._neuralNet = Pointer<NeuralNet>(new NeuralNet);
+            _allOrganisms[i]._neuralNet->create(cppn, *_nnPrms);
+        }
+    };
+
+    size_t threadCount = max(1u, Thread::hardware_concurrency());
+    size_t chunkSize   = _allOrganisms.size() / threadCount;
+
+    Vector<Thread> chunks(threadCount);
+
+    for (size_t i = 0; i < threadCount - 1; ++i) {
+        chunks[i] = Thread(generateChunk, i * chunkSize, (i + 1) * chunkSize);
+    }
+
+    chunks.back() = Thread(generateChunk, (chunks.size() - 1) * chunkSize, _allOrganisms.size());
+
+    for (auto& i : chunks) {
+        i.join();
+    }
+}
+
+void
+Population::replaceOrganism()
+{
+    _lastReplacement  = killPoorOrganism();
+    _lastMother       = nullptr;
+    _lastFather       = nullptr;
+    auto parentSpecie = chooseParentSpecie();
+    Vector<Organism*> eligibleMothers;
+
+    for (auto& i : *parentSpecie) {
+        if (i->_lifetime >= _prms._minimumLifetime) {
+            eligibleMothers.emplace_back(i);
+        }
+    }
+
+    size_t motherIdx = getRandSize(0, eligibleMothers.size() - 1);
+    _lastMother      = eligibleMothers[motherIdx];
+
+    if (getChance(_prms._sexualReproductionRate)) {
+        if (getChance(_prms._interspeciesMatingRate)) {
+            Vector<Organism*> eligibleFathers;
+
+            for (auto& i : _species) {
+                if (&i == parentSpecie) {
+                    continue;
+                }
+
+                for (auto& j : i) {
+                    if (j->_lifetime >= _prms._minimumLifetime) {
+                        eligibleFathers.emplace_back(j);
+                    }
+                }
+            }
+
+            if (eligibleFathers.size() > 0) {
+                size_t fatherIdx = getRandSize(0, eligibleFathers.size() - 1);
+                _lastFather      = eligibleFathers[fatherIdx];
+            }
+        } else {
+            if (eligibleMothers.size() > 1) {
+                size_t fatherIdx = getRandSize(0, eligibleMothers.size() - 2);
+
+                if (fatherIdx >= motherIdx) {
+                    ++fatherIdx;
+                }
+
+                _lastFather = eligibleMothers[fatherIdx];
+            }
+        }
+
+        if (!_lastFather) {
+            breedAsexually(_lastReplacement->_genome, _lastMother->_genome);
+        } else {
+            if (_lastFather->_fitness > _lastMother->_fitness) {
+                swap(_lastFather, _lastMother);
+            }
+
+            breedSexually(_lastReplacement->_genome, _lastMother->_genome, _lastFather->_genome);
+        }
+    } else {
+        breedAsexually(_lastReplacement->_genome, _lastMother->_genome);
+    }
+
+    assignToSpecie(*_lastReplacement);
+    ++_replacements;
+    --_oldOrganisms;
+
+    if ((_replacements % _prms._replBeforeReorganization) == 0) {
+        organizeSpecies();
+    }
+}
+
+Organism*
+Population::killPoorOrganism()
+{
+    Organism* poorOrganism  = nullptr;
+    double    minAdjFitness = -1.0;
+
+    for (auto& i : _allOrganisms) {
+        if (i._lifetime >= _prms._minimumLifetime && !i._isLocked) {
+            double orgAdjFitness = i._fitness / static_cast<double>(_species[i._specie].size());
+
+            if (orgAdjFitness < minAdjFitness || minAdjFitness == -1.0) {
+                minAdjFitness = orgAdjFitness;
+                poorOrganism  = &i;
+            }
+        }
+    }
+
+    if (poorOrganism->_isFrozen) {
+        poorOrganism->_isFrozen = false;
+        --_frozenOrganisms;
+    }
+
+    poorOrganism->_fitness  = 0.0;
+    poorOrganism->_lifetime = 0;
+    poorOrganism->_genome._nodeGenes.clear();
+
+    if (poorOrganism->_behavior) {
+        poorOrganism->_behavior->reset();
+    }
+
+    auto& childSpec = _species[poorOrganism->_specie];
+    auto  childIdx  = find(childSpec.begin(), childSpec.end(), poorOrganism);
+    childSpec.erase(childIdx);
+
+    return poorOrganism;
+}
+
+Vector<Organism*>*
+Population::chooseParentSpecie()
+{
+    double totalAverageFitnesses = 0.0;
+    Vector<double> averageFitnesses;
+    averageFitnesses.reserve(_species.size());
+
+    for (auto &i : _species) {
+        double specieFitness = 0.0;
+        size_t specieSize    = 0;
+
+        for (auto &j : i) {
+            if (j->_lifetime >= _prms._minimumLifetime) {
+                specieFitness += j->_fitness;
+                ++specieSize;
+            }
+        }
+
+        if (specieSize == 0) {
+            averageFitnesses.emplace_back(0.0);
+        } else {
+            specieFitness /= static_cast<double>(specieSize);
+            averageFitnesses.emplace_back(specieFitness);
+            totalAverageFitnesses += specieFitness;
+        }
+    }
+
+    size_t specieIdx = 0;
+    double selector  = getRandReal(0.0, totalAverageFitnesses) - averageFitnesses.front();
+
+    while (selector > 0.0) {
+        ++specieIdx;
+        selector -= averageFitnesses[specieIdx];
+    }
+
+    return &_species[specieIdx];
+}
+
+void
+Population::breedAsexually(Genome& child, const Genome& mother)
+{
+    child = mother;
+
+    if (mutateNodesAndLinks(child)) {
+        return;
+    }
+
+    for (auto& i : child._nodeGenes) {
+        for (auto& j : i.second._linkGenes) {
+            if (getChance(_prms._weightMutationRate)) {
+                j.second._weight += getWeightDeviation();
+
+                if (j.second._weight > _prms._weightRange) {
+                    j.second._weight = _prms._weightRange;
+                } else if (j.second._weight < -_prms._weightRange) {
+                    j.second._weight = -_prms._weightRange;
+                }
+            }
+        }
+    }
+}
+
+void
+Population::breedSexually(Genome& child, const Genome& mother, const Genome& father)
+{
+    child = mother;
+
+    for (auto& i : child._nodeGenes) {
+        if (father._nodeGenes.count(i.first)) {
+            auto& fatherNode = father._nodeGenes.at(i.first);
+
+            for (auto& j : i.second._linkGenes) {
+                if (fatherNode._linkGenes.count(j.first)) {
+                    auto& childLink  = j.second;
+                    auto& fatherLink = fatherNode._linkGenes.at(j.first);
+
+                    childLink._weight += fatherLink._weight;
+                    childLink._weight /= 2.0;
+
+                    if (!childLink._isEnabled || !fatherLink._isEnabled) {
+                        if (getChance(_prms._geneDisablingRatio)) {
+                            childLink._isEnabled = false;
+                        } else {
+                            childLink._isEnabled = true;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+bool
+Population::mutateNodesAndLinks(Genome& child)
+{
+    class LinkMutation {
+    public:
+        LinkMutation() = default;
+        LinkMutation(size_t source, size_t target)
+            : _source(source), _target(target)
+        {}
+
+        size_t _source = 0;
+        size_t _target = 0;
+    };
+
+    class NodeMutation {
+    public:
+        NodeMutation() = default;
+        NodeMutation(size_t source, size_t target, double weight)
+            : _source(source), _target(target), _weight(weight)
+        {}
+
+        size_t _source = 0;
+        size_t _target = 0;
+        double _weight = 0.0;
+    };
+
+    if (getChance(_prms._linkMutationRate)) {
+        Vector<LinkMutation> linkMutations;
+
+        for (auto& i : child._nodeGenes) {
+            auto& childNode  = i.first;
+            auto& childLinks = i.second._linkGenes;
+            auto& childDepth = i.second._depth;
+
+            for (size_t j = 0; j < child._inputs; ++j) {
+                if (!childLinks.count(j)) {
+                    linkMutations.emplace_back(j, childNode);
+                }
+            }
+
+            for (auto& j : child._nodeGenes) {
+                auto& sourceNode  = j.first;
+                auto& sourceDepth = j.second._depth;
+
+                if (!childLinks.count(sourceNode) && childDepth > sourceDepth) {
+                    linkMutations.emplace_back(sourceNode, childNode);
+                }
+            }
+        }
+
+        if (!linkMutations.empty()) {
+            size_t mutationIdx = getRandSize(0, linkMutations.size() - 1);
+            size_t target      = linkMutations[mutationIdx]._target;
+            size_t source      = linkMutations[mutationIdx]._source;
+
+            child._nodeGenes[target]._linkGenes[source] = {0.0};
+        }
+
+        return true;
+    }
+
+    if (getChance(_prms._nodeMutationRate)) {
+        Vector<NodeMutation> nodeMutations;
+
+        for (auto& i : child._nodeGenes) {
+            for (auto& j : i.second._linkGenes) {
+                auto& childLink = j.second;
+
+                if (childLink._isEnabled) {
+                    nodeMutations.emplace_back(j.first, i.first, childLink._weight);
+                }
+            }
+        }
+
+        size_t mutationIdx = getRandSize(0, nodeMutations.size() - 1);
+        size_t target      = nodeMutations[mutationIdx]._target;
+        size_t source      = nodeMutations[mutationIdx]._source;
+        double targetDepth = child._nodeGenes[target]._depth;
+        double sourceDepth = 0.0;
+
+        if (source >= child._inputs) {
+            sourceDepth = child._nodeGenes[source]._depth;
+        }
+
+        double   depth    = (targetDepth + sourceDepth) / 2.0;
+        NodeType function = getRandNodeType();
+
+        Innovation innov = {0, source, target, depth, function};
+        auto       iter  = find(_innovations.begin(), _innovations.end(), innov);
+        size_t     newID = 0;
+
+        if (iter != _innovations.end()) {
+            newID = iter->_number;
+        } else {
+            newID = _innovations.size() + _basicInnovs;
+            _innovations.emplace_back(newID, source, target, depth, function);
+        }
+
+        child._nodeGenes[newID]                                = {depth, function};
+        child._nodeGenes[newID]._linkGenes[source]             = {0.0};
+        child._nodeGenes[target]._linkGenes[newID]             = {0.0};
+        child._nodeGenes[target]._linkGenes[source]._isEnabled = false;
+
+        return true;
+    }
+
+    return false;
+}
+
+void
+Population::assignToSpecie(Organism& org)
+{
+    bool included = false;
+
+    for (auto i = _species.begin(), end = _species.end(); i != end; ++i) {
+        if (i->empty()) {
+            continue;
+        }
+
+        double genDistance = computeDistance((*i)[0]->_genome, org._genome);
+
+        if (genDistance < _distanceThreshold) {
+            included    = true;
+            org._specie = distance(_species.begin(), i);
+            i->emplace_back(&org);
+            break;
+        }
+    }
+
+    if (!included) {
+        org._specie = _species.size();
+        _species.emplace_back(1, &org);
+    }
+}
+
+void
+Population::organizeSpecies()
+{
+    _species.clear();
+
+    for (auto& i : _allOrganisms) {
+        assignToSpecie(i);
+    }
+
+    if (_species.size() < _prms._targetSpeciesCount && _distanceThreshold > _prms._distanceThresholdShift) {
+        _distanceThreshold -= _prms._distanceThresholdShift;
+    } else if (_species.size() > _prms._targetSpeciesCount) {
+        _distanceThreshold += _prms._distanceThresholdShift;
+    }
+}
+
+
+double
+Population::computeDistance(const Genome& g1, const Genome& g2) const
+{
+    class LinkPair {
+    public:
+        const Genome::NodeGene::LinkGene* _l1 = nullptr;
+        const Genome::NodeGene::LinkGene* _l2 = nullptr;
+    };
+
+    using LinkPairs = Map<size_t, LinkPair>;
+
+    class GenePair {
+    public:
+        const Genome::NodeGene* _g1 = nullptr;
+        const Genome::NodeGene* _g2 = nullptr;
+        LinkPairs _links;
+    };
+
+    using GenePairs = Map<size_t, GenePair>;
+
+    auto alignGenes = [](GenePairs& gPairs, const Genome& g1, const Genome& g2) {
+        for (auto& i : g1._nodeGenes) {
+            gPairs[i.first]._g1 = &i.second;
+
+            for (auto& j : i.second._linkGenes) {
+                gPairs[i.first]._links[j.first]._l1 = &j.second;
+            }
+        }
+
+        for (auto& i : g2._nodeGenes) {
+            gPairs[i.first]._g2 = &i.second;
+
+            for (auto& j : i.second._linkGenes) {
+                gPairs[i.first]._links[j.first]._l2 = &j.second;
+            }
+        }
+    };
+
+    GenePairs gPairs;
+    alignGenes(gPairs, g1, g2);
+
+    double g1Size           = 0.0;
+    double g2Size           = 0.0;
+    double disjointGenes    = 0.0;
+    double linkGenePairs    = 0.0;
+    double weightDifference = 0.0;
+
+    for (auto& i : gPairs) {
+        if (i.second._g1 && i.second._g2) {
+            ++g1Size;
+            ++g2Size;
+
+            for (auto& j : i.second._links) {
+                if (j.second._l1 && j.second._l2) {
+                    ++g1Size;
+                    ++g2Size;
+
+                    ++linkGenePairs;
+                    weightDifference += fabs(j.second._l1->_weight - j.second._l2->_weight);
+                } else if (j.second._l1) {
+                    ++g1Size;
+                    ++disjointGenes;
+                } else {
+                    ++g2Size;
+                    ++disjointGenes;
+                }
+            }
+        } else if (i.second._g1) {
+            double geneSize = static_cast<double>(1 + i.second._links.size());
+            g1Size         += geneSize;
+            disjointGenes  += geneSize;
+        } else {
+            double geneSize = static_cast<double>(1 + i.second._links.size());
+            g2Size         += geneSize;
+            disjointGenes  += geneSize;
+        }
+    }
+
+    double normalize      = g1Size > g2Size ? g1Size : g2Size;
+    double disjointFactor = (_prms._c1Disjoint * disjointGenes) / normalize;
+    double weightFactor   = 0.0;
+
+    if (linkGenePairs != 0.0) {
+        weightFactor  = (_prms._c3WeightDifference * weightDifference) / linkGenePairs;
+    }
+
+    return disjointFactor + weightFactor;
+}
+
+size_t
+Population::getRandSeed() const
+{
+    return chrono::system_clock::now().time_since_epoch().count();
+}
+
+size_t
+Population::getRandSize(size_t low, size_t hi)
+{
+    return IntDist(low, hi)(_randGen);
+}
+
+double
+Population::getRandReal(double low, double hi)
+{
+    return RealDist(low, hi)(_randGen);
+}
+
+double
+Population::getRandWeight()
+{
+    return _weightSelector(_randGen);
+}
+
+double
+Population::getWeightDeviation()
+{
+    return _weightDeviator(_randGen);
+}
+
+NodeType
+Population::getRandNodeType()
+{
+    return static_cast<NodeType>(_nodeTypeSelector(_randGen));
+}
+
+bool
+Population::getChance(double ratio)
+{
+    return _chanceSelector(_randGen) < ratio;
+}
diff --git a/src/QuadTree.cpp b/src/QuadTree.cpp
new file mode 100644
index 0000000..d6686ff
--- /dev/null
+++ b/src/QuadTree.cpp
@@ -0,0 +1,55 @@
+#include <HyperNeat/QuadTree.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+QuadTree::QuadTree(double segment, double x, double y)
+    : _segment(segment), _x(x), _y(y)
+{}
+
+double
+QuadTree::getSegment() const
+{
+    return _segment;
+}
+
+double
+QuadTree::getX() const
+{
+    return _x;
+}
+
+double
+QuadTree::getY() const
+{
+    return _y;
+}
+
+void
+QuadTree::subdivide(Function<bool(QuadTree*)> subdivider)
+{
+    if (subdivider(this)) {
+        double newSeg = _segment / 2.0;
+        _children.resize(4);
+        _children[0] = {newSeg, _x - newSeg, _y - newSeg};
+        _children[1] = {newSeg, _x + newSeg, _y - newSeg};
+        _children[2] = {newSeg, _x - newSeg, _y + newSeg};
+        _children[3] = {newSeg, _x + newSeg, _y + newSeg};
+
+        for (auto &i : _children) {
+            i.subdivide(subdivider);
+        }
+    }
+}
+
+void
+QuadTree::traverse(Function<void(const QuadTree*)> traverser) const
+{
+    if (!_children.empty()) {
+        for (auto &i : _children) {
+            i.traverse(traverser);
+        }
+    } else {
+        traverser(this);
+    }
+}
diff --git a/src/Utils/LoadFile.cpp b/src/Utils/LoadFile.cpp
new file mode 100644
index 0000000..3849f76
--- /dev/null
+++ b/src/Utils/LoadFile.cpp
@@ -0,0 +1,280 @@
+#include <HyperNeat/Population.hpp>
+#include <HyperNeat/NoveltyMetric.hpp>
+#include <HyperNeat/NeuralNetPrms.hpp>
+#include <HyperNeat/Utils/LoadFile.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+LoadFile::LoadFile(Istream& stream)
+    : _stream(stream)
+{
+    _stream.setf(ios::boolalpha);
+}
+
+void
+LoadFile::loadPopulation(Population& population)
+{
+    size_t  speciesCnt       = 0;
+    size_t  innovationsCount = 0;
+    ssize_t lastReplacement  = 0;
+    ssize_t lastMother       = 0;
+    ssize_t lastFather       = 0;
+    bool    hasNoveltyMetric = false;
+    bool    hasNeuralNets    = false;
+    auto&   prms             = population._prms;
+
+    nextPrm() >> population._populationLock;
+    nextPrm() >> population._lockedOrganisms;
+    nextPrm() >> population._frozenOrganisms;
+    nextPrm() >> innovationsCount;
+    nextPrm() >> population._basicInnovs;
+    nextPrm() >> lastReplacement;
+    nextPrm() >> lastMother;
+    nextPrm() >> lastFather;
+    nextPrm() >> population._replacements;
+    nextPrm() >> population._distanceThreshold;
+    nextPrm() >> population._oldOrganisms;
+    nextPrm() >> population._minOldOrganisms;
+    nextPrm() >> hasNoveltyMetric;
+    nextPrm() >> population._updates;
+    nextPrm() >> speciesCnt;
+    nextPrm() >> hasNeuralNets;
+
+    loadPopulationPrms(prms);
+
+    if (hasNeuralNets) {
+        population._nnPrms = Pointer<NeuralNetPrms>(new NeuralNetPrms);
+        loadNeuralNetPrms(*population._nnPrms);
+    }
+
+    population._innovations.resize(innovationsCount);
+
+    for (auto& i : population._innovations) {
+        String nodeStr;
+
+        nextPrm() >> i._number;
+        nextPrm() >> i._source;
+        nextPrm() >> i._target;
+        nextPrm() >> i._depth;
+        nextPrm() >> nodeStr;
+
+        i._nodeType = stringToNode(nodeStr);
+    }
+
+    population._allOrganisms.resize(prms._popSize, Organism(&population));
+    population._species.resize(speciesCnt);
+
+    for (auto& i : population._allOrganisms) {
+        loadOrganism(i);
+    }
+
+    population._lastReplacement = (lastReplacement == -1 ? nullptr: &population._allOrganisms[lastReplacement]);
+    population._lastMother      = (lastMother      == -1 ? nullptr: &population._allOrganisms[lastMother]);
+    population._lastFather      = (lastFather      == -1 ? nullptr: &population._allOrganisms[lastFather]);
+
+    for (auto& i : population._species) {
+        size_t specieSize = 0;
+        nextPrm() >> specieSize;
+
+        while (specieSize--) {
+            size_t organismIdx = 0;
+            nextArrayValue() >> organismIdx;
+
+            i.emplace_back(&population._allOrganisms[organismIdx]);
+        }
+    }
+
+    if (hasNoveltyMetric) {
+        NoveltyMetricPrms nmPrms;
+        loadNoveltyMetricPrms(nmPrms);
+
+        population.setNoveltyMetric(nmPrms);
+        loadNoveltyMetric(*population._noveltyMetric);
+    }
+
+    if (hasNeuralNets) {
+        population.generateAllNeuralNets();
+    }
+
+    if (prms._seed != 0) {
+        population._randGen.seed(prms._seed);
+    } else {
+        population._randGen.seed(population.getRandSeed());
+    }
+
+    population._weightDeviator = BellDist(0.0, prms._weightDeviation);
+    population._weightSelector = RealDist(-prms._weightRange, prms._weightRange);
+
+    population._organismsBeingGenerated = 0;
+}
+
+void
+LoadFile::loadPopulationPrms(PopulationPrms& prms)
+{
+    nextPrm() >> prms._popSize;
+    nextPrm() >> prms._cppnInputs;
+    nextPrm() >> prms._cppnOutputs;
+    nextPrm() >> prms._seed;
+    nextPrm() >> prms._weightRange;
+    nextPrm() >> prms._c1Disjoint;
+    nextPrm() >> prms._c3WeightDifference;
+    nextPrm() >> prms._initialDistanceThreshold;
+    nextPrm() >> prms._distanceThresholdShift;
+    nextPrm() >> prms._sexualReproductionRate;
+    nextPrm() >> prms._weightMutationRate;
+    nextPrm() >> prms._weightDeviation;
+    nextPrm() >> prms._interspeciesMatingRate;
+    nextPrm() >> prms._geneDisablingRatio;
+    nextPrm() >> prms._linkMutationRate;
+    nextPrm() >> prms._nodeMutationRate;
+    nextPrm() >> prms._targetSpeciesCount;
+    nextPrm() >> prms._eligibilityRatio;
+    nextPrm() >> prms._minimumLifetime;
+    nextPrm() >> prms._replBeforeReorganization;
+}
+
+void
+LoadFile::loadNeuralNetPrms(NeuralNetPrms& prms)
+{
+    size_t inputs  = 0;
+    size_t outputs = 0;
+
+    nextPrm() >> prms._testGridLevel;
+    nextPrm() >> prms._maxQuadTreeLevel;
+    nextPrm() >> prms._minQuadTreeLevel;
+    nextPrm() >> prms._bandPruningThreshold;
+    nextPrm() >> prms._varianceThreshold;
+    nextPrm() >> prms._divisionThreshold;
+    nextPrm() >> prms._iterations;
+    nextPrm() >> inputs;
+    nextPrm() >> outputs;
+
+    prms._inputMap.resize(inputs);
+    prms._outputMap.resize(outputs);
+
+    for (auto& i : prms._inputMap) {
+        nextArrayValue() >> i._x;
+        nextArrayValue() >> i._y;
+    }
+
+    for (auto& i : prms._outputMap) {
+        nextArrayValue() >> i._x;
+        nextArrayValue() >> i._y;
+    }
+}
+
+void
+LoadFile::loadNoveltyMetric(NoveltyMetric& noveltyMetric)
+{
+    for (auto& i : noveltyMetric._behaviors) {
+        nextPrm() >> i._criteriaReached;
+        nextPrm() >> i._noveltyScore;
+        nextPrm() >> i._toBeArchived;
+
+        for (auto& j : i._characterization) {
+            nextArrayValue() >> j;
+        }
+    }
+
+    size_t archiveSize = 0;
+    nextPrm() >> archiveSize;
+    noveltyMetric._archive.resize(archiveSize);
+
+    for (auto& i : noveltyMetric._archive) {
+        i.resize(noveltyMetric._prms._characterizationSize, 0.0);
+
+        for (auto& j : i) {
+            nextArrayValue() >> j;
+        }
+    }
+}
+
+void
+LoadFile::loadNoveltyMetricPrms(NoveltyMetricPrms& noveltyMetricPrms)
+{
+    nextPrm() >> noveltyMetricPrms._noveltyThreshold;
+    nextPrm() >> noveltyMetricPrms._referenceOrganisms;
+    nextPrm() >> noveltyMetricPrms._characterizationSize;
+    nextPrm() >> noveltyMetricPrms._criteriaReachedByDefault;
+}
+
+void
+LoadFile::loadOrganism(Organism& organism)
+{
+    nextPrm() >> organism._index;
+    nextPrm() >> organism._fitness;
+    nextPrm() >> organism._isLocked;
+    nextPrm() >> organism._isFrozen;
+    nextPrm() >> organism._specie;
+    nextPrm() >> organism._lifetime;
+
+    loadGenome(organism._genome);
+}
+
+void
+LoadFile::loadGenome(Genome& genome)
+{
+    size_t nodes = 0;
+
+    nextPrm() >> genome._inputs;
+    nextPrm() >> nodes;
+
+    while (nodes--) {
+        String nodeType;
+        size_t links = 0;
+        size_t innov = 0;
+
+        nextPrm() >> innov;
+        auto& nodeGene = genome._nodeGenes[innov];
+
+        nextPrm() >> nodeGene._depth;
+        nextPrm() >> nodeType;
+        nodeGene._nodeType = stringToNode(nodeType);
+
+        nextPrm() >> links;
+
+        while (links--) {
+            size_t source = 0;
+
+            nextPrm() >> source;
+            auto& linkGene = nodeGene._linkGenes[source];
+
+            nextPrm() >> linkGene._weight;
+            nextPrm() >> linkGene._isEnabled;
+        }
+    }
+}
+
+Istream&
+LoadFile::nextPrm(bool arrayVal)
+{
+    for (;;) {
+        char ch = static_cast<char>(_stream.peek());
+
+        if (ch == '#') {
+            _stream.ignore(numeric_limits<streamsize>::max(), '\n');
+        } else {
+            if (!arrayVal) {
+                _stream.ignore();
+
+                if (ch == '=') {
+                    return _stream;
+                }
+            } else {
+                if (string("-0123456789").find(ch) != string::npos) {
+                    return _stream;
+                }
+
+                _stream.ignore();
+            }
+        }
+    }
+}
+
+
+Istream&
+LoadFile::nextArrayValue()
+{
+    return nextPrm(true);
+}
diff --git a/src/Utils/NodeTypes.cpp b/src/Utils/NodeTypes.cpp
new file mode 100644
index 0000000..b11b079
--- /dev/null
+++ b/src/Utils/NodeTypes.cpp
@@ -0,0 +1,41 @@
+#include <HyperNeat/Utils/NodeTypes.hpp>
+
+namespace hyperneat
+{
+    String
+    nodeToString(NodeType type)
+    {
+        switch (type) {
+            case NodeType::SIGMOID:
+                return "\"sigmoid\"";
+
+            case NodeType::GAUSSIAN:
+                return "\"gaussian\"";
+
+            case NodeType::SINE:
+                return "\"sine\"";
+
+            case NodeType::ABSOLUTE:
+                return "\"absolute\"";
+
+            default:
+                return "\"nullType\"";
+        }
+    }
+
+    NodeType
+    stringToNode(const String& str)
+    {
+        if (str == "\"sigmoid\"") {
+            return NodeType::SIGMOID;
+        } else if (str == "\"gaussian\"") {
+            return NodeType::GAUSSIAN;
+        } else if (str == "\"sine\"") {
+            return NodeType::SINE;
+        } else if (str == "\"absolute\"") {
+            return NodeType::ABSOLUTE;
+        } else {
+            return NodeType::NULL_TYPE;
+        }
+    }
+}
diff --git a/src/Utils/Point.cpp b/src/Utils/Point.cpp
new file mode 100644
index 0000000..7dd9610
--- /dev/null
+++ b/src/Utils/Point.cpp
@@ -0,0 +1,33 @@
+#include <cmath>
+#include <HyperNeat/Utils/Point.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+Point::Point(double x, double y)
+    : _x(x), _y(y)
+{}
+
+double
+Point::distance(const Point& other) const
+{
+    double x = _x - other._x;
+    double y = _y - other._y;
+
+    return sqrt(x * x + y * y);
+}
+
+bool
+Point::operator==(const Point& other) const{
+    return (_x == other._x) && (_y == other._y);
+}
+
+bool
+Point::operator<(const Point& other) const
+{
+    if (_x == other._x) {
+        return _y < other._y;
+    } else {
+        return _x < other._x;
+    }
+}
diff --git a/src/Utils/SaveFile.cpp b/src/Utils/SaveFile.cpp
new file mode 100644
index 0000000..8c471f6
--- /dev/null
+++ b/src/Utils/SaveFile.cpp
@@ -0,0 +1,249 @@
+#include <HyperNeat/Population.hpp>
+#include <HyperNeat/NoveltyMetric.hpp>
+#include <HyperNeat/NeuralNetPrms.hpp>
+#include <HyperNeat/Utils/SaveFile.hpp>
+
+using namespace std;
+using namespace hyperneat;
+
+SaveFile::SaveFile(Ostream& stream)
+    : _stream(stream)
+{
+    _stream.setf(ios::boolalpha);
+    _stream.precision(numeric_limits<double>::digits10);
+}
+
+void
+SaveFile::savePopulation(Population& population, bool shuttedDown, size_t tabs, const String& prefix)
+{
+    print(tabs) << "# hyperneat::Population data:"                                                           << newl();
+    print(tabs) << "# -------------------------------------------------------------------------------------" << newl();
+    print(tabs) << "# DO NOT erase or change the order of any entry! Data is read back ENTRY by ENTRY."      << newl(2);
+
+    ssize_t lastReplacement = (population.getLastReplacement() ? population.getLastReplacement()->getIndex() : -1);
+    ssize_t lastMother      = (population.getLastMother() ? population.getLastMother()->getIndex() : -1);
+    ssize_t lastFather      = (population.getLastFather() ? population.getLastFather()->getIndex() : -1);
+
+    print(tabs) << "[" << prefix << "population]"                                                    << newl();
+    print(tabs + 1) << "populationLock    = " << (shuttedDown ? false : population.isLocked())       << newl();
+    print(tabs + 1) << "lockedOrganisms   = " << (shuttedDown ? 0 : population.getLockedOrganisms()) << newl();
+    print(tabs + 1) << "frozenOrganisms   = " << (shuttedDown ? 0 : population.getFrozenOrganisms()) << newl();
+    print(tabs + 1) << "innovationsCount  = " << population.getInnovationsCount()                    << newl();
+    print(tabs + 1) << "basicInnovations  = " << population.getBasicInnovationsCount()               << newl();
+    print(tabs + 1) << "lastReplacement   = " << lastReplacement                                     << newl();
+    print(tabs + 1) << "lastMother        = " << lastMother                                          << newl();
+    print(tabs + 1) << "lastFather        = " << lastFather                                          << newl();
+    print(tabs + 1) << "replacements      = " << population.getReplacements()                        << newl();
+    print(tabs + 1) << "distanceThreshold = " << population.getDistanceThreshold()                   << newl();
+    print(tabs + 1) << "oldOrganisms      = " << (shuttedDown ? 0 : population.getOldOrganisms())    << newl();
+    print(tabs + 1) << "minOldOrganisms   = " << population.getMinimumOldOrganisms()                 << newl();
+    print(tabs + 1) << "noveltyMetric     = " << population.isNoveltyMetricSet()                     << newl();
+    print(tabs + 1) << "updates           = " << population.getUpdates()                             << newl();
+    print(tabs + 1) << "species           = " << population.getSpecies().size()                      << newl();
+    print(tabs + 1) << "withNeuralNets    = " << population.hasNeuralNets()                          << newl(2);
+
+    savePopulationPrms(population.getPopulationPrms(), tabs + 1, prefix + "population.");
+
+    if (population.hasNeuralNets()) {
+        print() << newl();
+        saveNeuralNetPrms(population.getNeuralNetPrms(), tabs + 1, prefix + "population.");
+    }
+
+    for (auto& i : population.getInnovations()) {
+        print()                                                       << newl();
+        print(tabs + 1) << "[[" + prefix + "population.innovation]]"  << newl();
+        print(tabs + 2) << "number   = " << i._number                 << newl();
+        print(tabs + 2) << "source   = " << i._source                 << newl();
+        print(tabs + 2) << "target   = " << i._target                 << newl();
+        print(tabs + 2) << "depth    = " << i._depth                  << newl();
+        print(tabs + 2) << "function = " << nodeToString(i._nodeType) << newl();
+    }
+
+    for (auto& i : population.getAllOrganisms()) {
+        print() << newl();
+        saveOrganism(i, shuttedDown, tabs + 1, prefix + "population.");
+    }
+
+    for (auto& i : population.getSpecies()) {
+        print() << newl();
+        print(tabs + 1) << "[[" + prefix + "population.specie]]" << newl();
+        print(tabs + 2) << "size    = " << i.size()              << newl();
+        print(tabs + 2) << "members = ["                         << newl();
+
+        for (auto& j : i) {
+            print(tabs + 3) << j->getIndex() << "," << newl();
+        }
+
+        print(tabs + 2) << "]" << newl();
+    }
+
+    if (population.isNoveltyMetricSet()) {
+        print() << newl();
+        saveNoveltyMetric(population.getNoveltyMetric(), shuttedDown, tabs + 1, prefix + "population.");
+    }
+}
+
+void
+SaveFile::savePopulationPrms(const PopulationPrms& prms, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[" << prefix << "parameters]"                                      << newl();
+    print(tabs + 1) << "popSize                  = " << prms._popSize                  << newl();
+    print(tabs + 1) << "cppnInputs               = " << prms._cppnInputs               << newl();
+    print(tabs + 1) << "cppnOutputs              = " << prms._cppnOutputs              << newl();
+    print(tabs + 1) << "seed                     = " << prms._seed                     << newl();
+    print(tabs + 1) << "weightRange              = " << prms._weightRange              << newl();
+    print(tabs + 1) << "disjointCoeff            = " << prms._c1Disjoint               << newl();
+    print(tabs + 1) << "weightDifferenceCoeff    = " << prms._c3WeightDifference       << newl();
+    print(tabs + 1) << "initialDistanceThreshold = " << prms._initialDistanceThreshold << newl();
+    print(tabs + 1) << "distanceThresholdShift   = " << prms._distanceThresholdShift   << newl();
+    print(tabs + 1) << "sexualReproductionRate   = " << prms._sexualReproductionRate   << newl();
+    print(tabs + 1) << "weightMutationRate       = " << prms._weightMutationRate       << newl();
+    print(tabs + 1) << "weightDeviation          = " << prms._weightDeviation          << newl();
+    print(tabs + 1) << "interspeciesMatingRate   = " << prms._interspeciesMatingRate   << newl();
+    print(tabs + 1) << "geneDisablingRatio       = " << prms._geneDisablingRatio       << newl();
+    print(tabs + 1) << "linkMutationRate         = " << prms._linkMutationRate         << newl();
+    print(tabs + 1) << "nodeMutationRate         = " << prms._nodeMutationRate         << newl();
+    print(tabs + 1) << "targetSpeciesCount       = " << prms._targetSpeciesCount       << newl();
+    print(tabs + 1) << "eligibilityRatio         = " << prms._eligibilityRatio         << newl();
+    print(tabs + 1) << "minimumLifetime          = " << prms._minimumLifetime          << newl();
+    print(tabs + 1) << "replBeforeReorganization = " << prms._replBeforeReorganization << newl();
+}
+
+void
+SaveFile::saveNeuralNetPrms(const NeuralNetPrms& prms, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[" << prefix << "neuralNetParameters]"                     << newl();
+    print(tabs + 1) << "testGridLevel        = " << prms._testGridLevel        << newl();
+    print(tabs + 1) << "maxQuadTreeLevel     = " << prms._maxQuadTreeLevel     << newl();
+    print(tabs + 1) << "minQuadTreeLevel     = " << prms._minQuadTreeLevel     << newl();
+    print(tabs + 1) << "bandPruningThreshold = " << prms._bandPruningThreshold << newl();
+    print(tabs + 1) << "varianceThreshold    = " << prms._varianceThreshold    << newl();
+    print(tabs + 1) << "divisionThreshold    = " << prms._divisionThreshold    << newl();
+    print(tabs + 1) << "iterations           = " << prms._iterations           << newl();
+    print(tabs + 1) << "inputs               = " << prms._inputMap.size()      << newl();
+    print(tabs + 1) << "outputs              = " << prms._outputMap.size()     << newl(2);
+    print(tabs + 1) << "inputMap = ["                                          << newl();
+
+    for (auto& i : prms._inputMap) {
+        print(tabs + 2) << "[" << i._x << ", " << i._y << "]," << newl();
+    }
+
+    print(tabs + 1) << "]"             << newl(2);
+    print(tabs + 1) << "outputMap = [" << newl();
+
+    for (auto& i : prms._outputMap) {
+        print(tabs + 2) << "[" << i._x << ", " << i._y << "]," << newl();
+    }
+
+    print(tabs + 1) << "]" << newl();
+}
+
+void
+SaveFile::saveOrganism(const Organism& organism, bool shuttedDown, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[[" << prefix << "organism]]" << newl();
+
+    if (organism.isChampion()) {
+        print(tabs + 1) << "# Current champion." << newl(2);
+    }
+
+    print(tabs + 1) << "index    = " << organism.getIndex()                         << newl();
+    print(tabs + 1) << "fitness  = " << (shuttedDown ? 0 : organism._fitness)       << newl();
+    print(tabs + 1) << "isLocked = " << (shuttedDown ? false : organism.isLocked()) << newl();
+    print(tabs + 1) << "isFrozen = " << (shuttedDown ? false : organism.isFrozen()) << newl();
+    print(tabs + 1) << "specie   = " << organism.getSpecie()                        << newl();
+    print(tabs + 1) << "lifetime = " << (shuttedDown ? 0 : organism.getLifetime())  << newl(2);
+
+    saveGenome(organism.getGenome(), tabs + 1, prefix + "organism.");
+
+    // If shutted down, save Neural Net
+}
+
+void
+SaveFile::saveGenome(const Genome& genome, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[" << prefix << "genome]"                     << newl();
+    print(tabs + 1) << "inputs = "    << genome._inputs           << newl();
+    print(tabs + 1) << "nodes  = "    << genome._nodeGenes.size() << newl();
+
+    for (auto& i : genome._nodeGenes) {
+        print()                                                                << newl();
+        print(tabs + 1) << "[[" << prefix << "genome.nodeGene]]"               << newl();
+        print(tabs + 2) << "innovation = " << i.first                          << newl();
+        print(tabs + 2) << "depth      = " << i.second._depth                  << newl();
+        print(tabs + 2) << "function   = " << nodeToString(i.second._nodeType) << newl();
+        print(tabs + 2) << "links      = " << i.second._linkGenes.size()       << newl();
+
+        for (auto& j : i.second._linkGenes) {
+            print()                                                           << newl();
+            print(tabs + 2) << "[[" << prefix << "genome.nodeGene.linkGene]]" << newl();
+            print(tabs + 3) << "source    = " << j.first                      << newl();
+            print(tabs + 3) << "weight    = " << j.second._weight             << newl();
+            print(tabs + 3) << "isEnabled = " << j.second._isEnabled          << newl();
+        }
+    }
+}
+
+void
+SaveFile::saveNoveltyMetric(const NoveltyMetric& noveltyMetric, bool shuttedDown, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[" << prefix << "noveltyMetric]" << newl();
+    saveNoveltyMetricPrms(noveltyMetric.getPrms(), tabs + 1, prefix + "noveltyMetric.");
+
+    bool critDef = noveltyMetric.getPrms()._criteriaReachedByDefault;
+
+    for (auto& i : noveltyMetric.getBehaviors()) {
+        print()                                                                                   << newl();
+        print(tabs + 1) << "[[" << prefix << "noveltyMetric.behavior]]"                           << newl();
+        print(tabs + 2) << "criteriaReached  = " << (shuttedDown ? critDef : i._criteriaReached)  << newl();
+        print(tabs + 2) << "noveltyScore     = " << (shuttedDown ? 0 : i.getNoveltyScore())       << newl();
+        print(tabs + 2) << "toBeArchived     = " << (shuttedDown ? false : i.isToBeArchived())    << newl();
+        print(tabs + 2) << "characterization = ["                                                 << newl();
+
+        for (auto& j : i.getCharacterization()) {
+            print(tabs + 3) << j << "," << newl();
+        }
+
+        print(tabs + 2) << "]" << newl();
+    }
+
+    print()                                                                        << newl();
+    print(tabs + 1) << "[" << prefix << "noveltyMetric.archive]"                   << newl();
+    print(tabs + 2) << "size              = " << noveltyMetric.getArchive().size() << newl();
+    print(tabs + 2) << "characterizations = ["                                     << newl();
+
+    for (auto& i : noveltyMetric.getArchive()) {
+        print(tabs + 3) << "[" << newl();
+
+        for (auto& j : i) {
+            print(tabs + 4) << j << "," << newl();
+        }
+
+        print(tabs + 3) << "]," << newl();
+    }
+
+    print(tabs + 2) << "]" << newl();
+}
+
+void
+SaveFile::saveNoveltyMetricPrms(const NoveltyMetricPrms& noveltyMetricPrms, size_t tabs, const String& prefix)
+{
+    print(tabs) << "[" << prefix << "parameters]" << newl();
+    print(tabs + 1) << "noveltyThreshold         = " << noveltyMetricPrms._noveltyThreshold         << newl();
+    print(tabs + 1) << "referenceOrganisms       = " << noveltyMetricPrms._referenceOrganisms       << newl();
+    print(tabs + 1) << "characterizationSize     = " << noveltyMetricPrms._characterizationSize     << newl();
+    print(tabs + 1) << "criteriaReachedByDefault = " << noveltyMetricPrms._criteriaReachedByDefault << newl();
+}
+
+Ostream&
+SaveFile::print(size_t tabs)
+{
+    _stream << String(tabs * 4, ' ');
+    return _stream;
+}
+
+String
+SaveFile::newl(size_t lines)
+{
+    return String(lines, '\n');
+}
diff --git a/src/Utils/ValuePoint.cpp b/src/Utils/ValuePoint.cpp
new file mode 100644
index 0000000..d9cdc64
--- /dev/null
+++ b/src/Utils/ValuePoint.cpp
@@ -0,0 +1,7 @@
+#include <HyperNeat/Utils/ValuePoint.hpp>
+
+using namespace hyperneat;
+
+ValuePoint::ValuePoint(double x, double y, double value, double segment)
+    : Point(x, y), _value(value), _segment(segment)
+{}