Hadrons: first version of the genetic scheduler

programs/Hadrons/Application.cc

@@ -26,7 +26,7 @@ directory.
 *******************************************************************************/
 
 #include <Hadrons/Application.hpp>
-#include <Hadrons/Graph.hpp>
+#include <Hadrons/GeneticScheduler.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -115,7 +115,7 @@ sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
 
 void Application::schedule(void)
 {
-    // memory peak and invers memory peak functions
+    // memory peak function
     auto memPeak = [this](const std::vector<std::string> &program)
     {
         unsigned int memPeak;
@@ -131,16 +131,11 @@ void Application::schedule(void)
         
         return memPeak;
     };
-    auto invMemPeak = [&memPeak](const std::vector<std::string> &program)
-    {
-        return 1./memPeak(program);
-    };
-    
-    Graph<std::string>            moduleGraph;
-    
-    LOG(Message) << "Scheduling computation..." << std::endl;
     
     // create dependency graph
+    Graph<std::string> moduleGraph;
+    
+    LOG(Message) << "Scheduling computation..." << std::endl;
     for (auto &m: module_)
     {
         std::vector<std::string> input = m.second->getInput();
@@ -157,39 +152,47 @@ void Application::schedule(void)
         }
     }
     
-    // topological sort
+    // constrained topological sort using a genetic algorithm
-    unsigned int k = 0;
+    constexpr unsigned int maxGen = 200, maxCstGen = 50;
-    
+    unsigned int           k = 0, gen, prevPeak, nCstPeak = 0;
     std::vector<Graph<std::string>> con = moduleGraph.getConnectedComponents();
-    
+    GeneticScheduler<std::string>::Parameters par;
+    std::random_device rd;
+
+    par.popSize      = 20;
+    par.mutationRate = .1;
+    par.seed         = rd();
+    CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed));
     for (unsigned int i = 0; i < con.size(); ++i)
     {
-//        std::vector<std::vector<std::string>> t = con[i].allTopoSort();
+        GeneticScheduler<std::string> scheduler(con[i], memPeak, par);
-//        int                                   memPeak, minMemPeak = -1;
+        
-//        unsigned int                          bestInd;
+        gen = 0;
-//        bool                                  msg;
+        do
-//
+        {
-//        LOG(Message) << "analyzing " << t.size() << " possible programs..."
+            scheduler.nextGeneration();
-//                     << std::endl;
+            if (gen != 0)
-//        env_.dryRun(true);
+            {
-//        for (unsigned int p = 0; p < t.size(); ++p)
+                if (prevPeak == scheduler.getMinValue())
-//        {
+                {
-//            msg = HadronsLogMessage.isActive();
+                    nCstPeak++;
-//            HadronsLogMessage.Active(false);
+                }
-//        
+                else
-//            memPeak = execute(t[p]);
+                {
-//            if ((memPeak < minMemPeak) or (minMemPeak < 0))
+                    nCstPeak = 0;
-//            {
+                }
-//                minMemPeak = memPeak;
+            }
-//                bestInd = p;
+            prevPeak = scheduler.getMinValue();
-//            }
+            if (gen % 10 == 0)
-//            HadronsLogMessage.Active(msg);
+            {
-//            env_.freeAll();
+                LOG(Iterative) << "Generation " << gen << ": "
-//        }
+                               << MEM_MSG(scheduler.getMinValue()) << std::endl;
-//        env_.dryRun(false);
+            }
-        std::vector<std::string> t = con[i].topoSort();
+            gen++;
+        } while ((gen < maxGen) and (nCstPeak < maxCstGen));
+        auto &t = scheduler.getMinSchedule();
         LOG(Message) << "Program " << i + 1 << " (memory peak: "
-                     << MEM_MSG(memPeak(t)) << "):" << std::endl;
+                     << MEM_MSG(scheduler.getMinValue()) << "):" << std::endl;
         for (unsigned int j = 0; j < t.size(); ++j)
         {
             program_.push_back(t[j]);
@@ -213,13 +216,14 @@ void Application::configLoop(void)
         execute(program_);
         env_.freeAll();
     }
+    LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
 }
 
 unsigned int Application::execute(const std::vector<std::string> &program)
 {
-    unsigned int                       memPeak = 0, size;
+    unsigned int                       memPeak = 0, sizeBefore, sizeAfter;
     std::vector<std::set<std::string>> freeProg;
-    bool                               continueCollect;
+    bool                               continueCollect, nothingFreed;
     
     // build garbage collection schedule
     freeProg.resize(program.size());
@@ -247,15 +251,17 @@ unsigned int Application::execute(const std::vector<std::string> &program)
                      << program.size() << " (module '" << program[i] << "') "
                      << SEP << std::endl;
         (*module_[program[i]])();
-        size = env_.getTotalSize();
+        sizeBefore = env_.getTotalSize();
         // print used memory after execution
-        LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl;
+        LOG(Message) << "Allocated objects: " << MEM_MSG(sizeBefore)
-        if (size > memPeak)
+                     << std::endl;
+        if (sizeBefore > memPeak)
         {
-            memPeak = size;
+            memPeak = sizeBefore;
         }
         // garbage collection for step i
         LOG(Message) << "Garbage collection..." << std::endl;
+        nothingFreed = true;
         do
         {
             continueCollect = false;
@@ -270,6 +276,7 @@ unsigned int Application::execute(const std::vector<std::string> &program)
                     // if an object has been freed, remove it from
                     // the garbage collection schedule
                     freeProg[i].erase(n);
+                    nothingFreed = false;
                 }
             }
         } while (continueCollect);
@@ -282,9 +289,17 @@ unsigned int Application::execute(const std::vector<std::string> &program)
                 freeProg[i + 1].insert(n);
             }
         }
-        // print used memory after garbage collection
+        // print used memory after garbage collection if necessary
-        size = env_.getTotalSize();
+        sizeAfter = env_.getTotalSize();
-        LOG(Message) << "Allocated objects: " << MEM_MSG(size) << std::endl;
+        if (sizeBefore != sizeAfter)
+        {
+            LOG(Message) << "Allocated objects: " << MEM_MSG(sizeAfter)
+                         << std::endl;
+        }
+        else
+        {
+            LOG(Message) << "Nothing to free" << std::endl;
+        }
     }
     
     return memPeak;

programs/Hadrons/CMeson.cc

@@ -64,7 +64,7 @@ std::vector<std::string> CMeson::getOutput(void)
 void CMeson::execute(void)
 {
     LOG(Message) << "Computing meson contraction '" << getName() << "' using"
-                 << " quarks '" << par_.q1 << " and '" << par_.q2 << "'"
+                 << " quarks '" << par_.q1 << "' and '" << par_.q2 << "'"
                  << std::endl;
     
     XmlWriter             writer(par_.output);

programs/Hadrons/Environment.cc

@@ -69,10 +69,13 @@ unsigned int Environment::getTrajectory(void) const
 // grids ///////////////////////////////////////////////////////////////////////
 void Environment::createGrid(const unsigned int Ls)
 {
-    auto g = getGrid();
+    if (grid5d_.find(Ls) == grid5d_.end())
-    
+    {
-    grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
+        auto g = getGrid();
-    gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
+        
+        grid5d_[Ls].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
+        gridRb5d_[Ls].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
+    }
 }
 
 GridCartesian * Environment::getGrid(const unsigned int Ls) const
@@ -329,11 +332,14 @@ void Environment::addOwnership(const std::string owner,
 
 bool Environment::hasOwners(const std::string name) const
 {
-    try
+    
+    auto it = owners_.find(name);
+    
+    if (it != owners_.end())
     {
-        return (!owners_.at(name).empty());
+        return (!it->second.empty());
     }
-    catch (std::out_of_range &)
+    else
     {
         return false;
     }

programs/Hadrons/GeneticScheduler.hpp

@@ -0,0 +1,233 @@
+/*******************************************************************************
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: programs/Hadrons/GeneticScheduler.hpp
+
+Copyright (C) 2016
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution 
+directory.
+*******************************************************************************/
+
+#ifndef Hadrons_GeneticScheduler_hpp_
+#define Hadrons_GeneticScheduler_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Graph.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                   Scheduler based on a genetic algorithm                   *
+ ******************************************************************************/
+template <typename T>
+class GeneticScheduler
+{
+public:
+    typedef std::function<int(const std::vector<T> &)> ObjFunc;
+    struct Parameters
+    {
+        double       mutationRate;
+        unsigned int popSize, seed;
+    };
+public:
+    // constructor
+    GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                     const Parameters &par);
+    // destructor
+    virtual ~GeneticScheduler(void) = default;
+    // access
+    const std::vector<T> & getMinSchedule(void);
+    int                    getMinValue(void);
+    // breed a new generation
+    void nextGeneration(void);
+    // print population
+    friend std::ostream & operator<<(std::ostream &out,
+                                     const GeneticScheduler<T> &s)
+    {
+        for (auto &p: s.population_)
+        {
+            out << p.second << ": " << p.first << std::endl;
+        }
+        
+        return out;
+    }
+private:
+    // randomly initialize population
+    void initPopulation(void);
+    // genetic operators
+    const std::vector<T> & selection(void);
+    void                   crossover(const std::vector<T> &c1,
+                                     const std::vector<T> &c2);
+    void                   mutation(std::vector<T> &c);
+private:
+    Graph<T>                           &graph_;
+    const ObjFunc                      &func_;
+    const Parameters                   par_;
+    std::multimap<int, std::vector<T>> population_;
+    std::mt19937                       gen_;
+};
+
+/******************************************************************************
+ *                       template implementation                              *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename T>
+GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+                                      const Parameters &par)
+: graph_(graph)
+, func_(func)
+, par_(par)
+{
+    gen_.seed(par_.seed);
+}
+
+// access //////////////////////////////////////////////////////////////////////
+template <typename T>
+const std::vector<T> & GeneticScheduler<T>::getMinSchedule(void)
+{
+    return population_.begin()->second;
+}
+
+template <typename T>
+int GeneticScheduler<T>::getMinValue(void)
+{
+    return population_.begin()->first;
+}
+
+// breed a new generation //////////////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::nextGeneration(void)
+{
+    std::uniform_real_distribution<double> dis(0., 1.);
+    
+    // random initialization of the population if necessary
+    if (population_.size() != par_.popSize)
+    {
+        initPopulation();
+    }
+
+    // mating
+    for (unsigned int i = 0; i < par_.popSize/2; ++i)
+    {
+        auto &p1 = selection(), &p2 = selection();
+        crossover(p1, p2);
+    }
+    
+    // random mutations
+    auto buf = population_;
+    population_.clear();
+    for (auto &c: buf)
+    {
+        if (dis(gen_) < par_.mutationRate)
+        {
+            mutation(c.second);
+        }
+        population_.emplace(func_(c.second), c.second);
+    }
+    
+    // grim reaper
+    auto it = population_.begin();
+    
+    std::advance(it, par_.popSize);
+    population_.erase(it, population_.end());
+}
+
+// randomly initialize population //////////////////////////////////////////////
+template <typename T>
+void GeneticScheduler<T>::initPopulation(void)
+{
+    population_.clear();
+    for (unsigned int i = 0; i < par_.popSize; ++i)
+    {
+        auto p = graph_.topoSort(gen_);
+        
+        population_.emplace(func_(p), p);
+    }
+}
+
+// genetic operators ///////////////////////////////////////////////////////////
+template <typename T>
+const std::vector<T> & GeneticScheduler<T>::selection(void)
+{
+    std::vector<double> prob;
+    
+    for (auto &c: population_)
+    {
+        prob.push_back(1./c.first);
+    }
+    std::discrete_distribution<unsigned int> dis(prob.begin(), prob.end());
+    auto rIt = population_.begin();
+    std::advance(rIt, dis(gen_));
+    
+    return rIt->second;
+}
+
+template <typename T>
+void GeneticScheduler<T>::crossover(const std::vector<T> &p1,
+                                    const std::vector<T> &p2)
+{
+    std::uniform_int_distribution<unsigned int> dis(1, p1.size() - 2);
+    unsigned int                                cut = dis(gen_);
+    std::vector<T>                              c1, c2, buf;
+    
+    auto cross = [&buf, cut](std::vector<T> &c, const std::vector<T> &p1,
+                             const std::vector<T> &p2)
+    {
+        buf = p2;
+        for (unsigned int i = 0; i < cut; ++i)
+        {
+            c.push_back(p1[i]);
+            buf.erase(std::find(buf.begin(), buf.end(), p1[i]));
+        }
+        for (unsigned int i = 0; i < buf.size(); ++i)
+        {
+            c.push_back(buf[i]);
+        }
+    };
+    
+    cross(c1, p1, p2);
+    cross(c2, p2, p1);
+    population_.emplace(func_(c1), c1);
+    population_.emplace(func_(c2), c2);
+}
+
+template <typename T>
+void GeneticScheduler<T>::mutation(std::vector<T> &c)
+{
+    std::uniform_int_distribution<unsigned int> dis(1, c.size() - 2);
+    unsigned int                                cut = dis(gen_);
+    Graph<T>                                    g = graph_;
+    std::vector<T>                              buf;
+    
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        g.removeVertex(c[i]);
+    }
+    buf = g.topoSort(gen_);
+    for (unsigned int i = cut; i < c.size(); ++i)
+    {
+        buf.push_back(c[i]);
+    }
+    c = buf;
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_GeneticScheduler_hpp_

programs/Hadrons/Global.cc

@@ -34,4 +34,5 @@ using namespace Hadrons;
 HadronsLogger Hadrons::HadronsLogError(1,"Error");
 HadronsLogger Hadrons::HadronsLogWarning(1,"Warning");
 HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
+HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
 HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");

programs/Hadrons/Global.hpp

@@ -63,6 +63,7 @@ abort();
 extern HadronsLogger HadronsLogError;
 extern HadronsLogger HadronsLogWarning;
 extern HadronsLogger HadronsLogMessage;
+extern HadronsLogger HadronsLogIterative;
 extern HadronsLogger HadronsLogDebug;
 
 // singleton pattern

programs/Hadrons/Graph.hpp

@@ -56,13 +56,14 @@ public:
     // destructor
     virtual ~Graph(void) = default;
     // access
-    void         addVertex(const T &value);
+    void           addVertex(const T &value);
-    void         addEdge(const Edge &e);
+    void           addEdge(const Edge &e);
-    void         addEdge(const T &start, const T &end);
+    void           addEdge(const T &start, const T &end);
-    void         removeVertex(const T &value);
+    std::vector<T> getVertices(void) const;
-    void         removeEdge(const Edge &e);
+    void           removeVertex(const T &value);
-    void         removeEdge(const T &start, const T &end);
+    void           removeEdge(const Edge &e);
-    unsigned int size(void) const;
+    void           removeEdge(const T &start, const T &end);
+    unsigned int   size(void) const;
     // tests
     bool gotValue(const T &value) const;
     // graph topological manipulations
@@ -71,7 +72,9 @@ public:
     std::vector<T>              getParents(const T &value) const;
     std::vector<T>              getRoots(void) const;
     std::vector<Graph<T>>       getConnectedComponents(void) const;
-    std::vector<T>              topoSort(const bool randomize = false);
+    std::vector<T>              topoSort(void);
+    template <typename Gen>
+    std::vector<T>              topoSort(Gen &gen);
     std::vector<std::vector<T>> allTopoSort(void);
     // I/O
     friend std::ostream & operator<<(std::ostream &out, const Graph<T> &g)
@@ -97,9 +100,11 @@ private:
     void      unmarkAll(void);
     bool      isMarked(const T &value) const;
     const T * getFirstMarked(const bool isMarked = true) const;
-    const T * getRandomMarked(const bool isMarked = true);
+    template <typename Gen>
+    const T * getRandomMarked(const bool isMarked, Gen &gen);
     const T * getFirstUnmarked(void) const;
-    const T * getRandomUnmarked(void);
+    template <typename Gen>
+    const T * getRandomUnmarked(Gen &gen);
     // prune marked/unmarked vertices
     void removeMarked(const bool isMarked = true);
     void removeUnmarked(void);
@@ -109,7 +114,6 @@ private:
 private:
     std::map<T, bool>  isMarked_;
     std::set<Edge>     edgeSet_;
-    std::mt19937       gen_;
 };
 
 // build depedency matrix from topological sorts
@@ -127,11 +131,7 @@ makeDependencyMatrix(const std::vector<std::vector<T>> &topSort);
 // constructor /////////////////////////////////////////////////////////////////
 template <typename T>
 Graph<T>::Graph(void)
-{
+{}
-    std::random_device rd;
-    
-    gen_.seed(rd());
-}
 
 // access //////////////////////////////////////////////////////////////////////
 // complexity: log(V)
@@ -157,6 +157,19 @@ void Graph<T>::addEdge(const T &start, const T &end)
     addEdge(Edge(start, end));
 }
 
+template <typename T>
+std::vector<T> Graph<T>::getVertices(void) const
+{
+    std::vector<T> vertex;
+    
+    for (auto &v: isMarked_)
+    {
+        vertex.push_back(v.first);
+    }
+    
+    return vertex;
+}
+
 // complexity: O(V*log(V))
 template <typename T>
 void Graph<T>::removeVertex(const T &value)
@@ -311,7 +324,8 @@ const T * Graph<T>::getFirstMarked(const bool isMarked) const
 
 // complexity: O(log(V))
 template <typename T>
-const T * Graph<T>::getRandomMarked(const bool isMarked)
+template <typename Gen>
+const T * Graph<T>::getRandomMarked(const bool isMarked, Gen &gen)
 {
     auto pred = [&isMarked](const std::pair<T, bool> &v)
     {
@@ -320,7 +334,7 @@ const T * Graph<T>::getRandomMarked(const bool isMarked)
     std::uniform_int_distribution<unsigned int> dis(0, size() - 1);
     auto                                        rIt = isMarked_.begin();
     
-    std::advance(rIt, dis(gen_));
+    std::advance(rIt, dis(gen));
     auto vIt = std::find_if(rIt, isMarked_.end(), pred);
     if (vIt != isMarked_.end())
     {
@@ -349,9 +363,10 @@ const T * Graph<T>::getFirstUnmarked(void) const
 
 // complexity: O(log(V))
 template <typename T>
-const T * Graph<T>::getRandomUnmarked(void)
+template <typename Gen>
+const T * Graph<T>::getRandomUnmarked(Gen &gen)
 {
-    return getRandomMarked(false);
+    return getRandomMarked(false, gen);
 }
 
 // prune marked/unmarked vertices //////////////////////////////////////////////
@@ -515,7 +530,7 @@ std::vector<Graph<T>> Graph<T>::getConnectedComponents(void) const
 // topological sort using a directed DFS algorithm
 // complexity: O(V*log(V))
 template <typename T>
-std::vector<T> Graph<T>::topoSort(const bool randomize)
+std::vector<T> Graph<T>::topoSort(void)
 {
     std::stack<T>     buf;
     std::vector<T>    res;
@@ -534,10 +549,6 @@ std::vector<T> Graph<T>::topoSort(const bool randomize)
             std::vector<T> child = getChildren(v);
 
             tmpMarked[v] = true;
-            if (randomize)
-            {
-                std::shuffle(child.begin(), child.end(), gen_);
-            }
             for (auto &c: child)
             {
                 visit(c);
@@ -556,25 +567,71 @@ std::vector<T> Graph<T>::topoSort(const bool randomize)
     
     // loop on unmarked vertices
     unmarkAll();
-    if (randomize)
+    vPt = getFirstUnmarked();
-    {
-        vPt = getRandomUnmarked();
-    }
-    else
-    {
-        vPt = getFirstUnmarked();
-    }
     while (vPt)
     {
         visit(*vPt);
-        if (randomize)
+        vPt = getFirstUnmarked();
+    }
+    unmarkAll();
+    
+    // create result vector
+    while (!buf.empty())
+    {
+        res.push_back(buf.top());
+        buf.pop();
+    }
+    
+    return res;
+}
+
+// random version of the topological sort
+// complexity: O(V*log(V))
+template <typename T>
+template <typename Gen>
+std::vector<T> Graph<T>::topoSort(Gen &gen)
+{
+    std::stack<T>     buf;
+    std::vector<T>    res;
+    const T           *vPt;
+    std::map<T, bool> tmpMarked(isMarked_);
+    
+    // visit function
+    std::function<void(const T &)> visit = [&](const T &v)
+    {
+        if (tmpMarked.at(v))
         {
-            vPt = getRandomUnmarked();
+            HADRON_ERROR("cannot topologically sort a cyclic graph");
         }
-        else
+        if (!isMarked(v))
         {
-            vPt = getFirstUnmarked();
+            std::vector<T> child = getChildren(v);
+            
+            tmpMarked[v] = true;
+            std::shuffle(child.begin(), child.end(), gen);
+            for (auto &c: child)
+            {
+                visit(c);
+            }
+            mark(v);
+            tmpMarked[v] = false;
+            buf.push(v);
         }
+    };
+    
+    // reset temporary marks
+    for (auto &v: tmpMarked)
+    {
+        tmpMarked.at(v.first) = false;
+    }
+    
+    // loop on unmarked vertices
+    unmarkAll();
+    vPt = getRandomUnmarked(gen);
+    while (vPt)
+    {
+        visit(*vPt);
+        vPt = getRandomUnmarked(gen);
     }
     unmarkAll();
     

programs/Hadrons/Hadrons.cc

@@ -49,6 +49,7 @@ int main(int argc, char *argv[])
     HadronsLogError.Active(GridLogError.isActive());
     HadronsLogWarning.Active(GridLogWarning.isActive());
     HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
     HadronsLogDebug.Active(GridLogDebug.isActive());
     LOG(Message) << "Grid initialized" << std::endl;