parallel.cpp

Go to the documentation of this file.

#include <parallel.h>
#include <concurrency.h>

using namespace frechet;
using namespace poly;
using namespace reach;
using namespace app;

//  Algorithm factory
Algorithm::ptr frechet::poly::newPolyAlgorithm(bool topo) 
{    
    if (topo)
        return poly::Algorithm::ptr(new poly::AlgorithmTopoSort());
    else if (ConcurrencyContext::countThreads() <= 1)
        return poly::Algorithm::ptr(new poly::AlgorithmSingleCore());
    else
        return poly::Algorithm::ptr(new poly::AlgorithmMultiCore());
}

AlgorithmSingleCore::AlgorithmSingleCore()
:   Algorithm()
{ }

AlgorithmMultiCore::LocalTriangulation AlgorithmMultiCore::qtri;
AlgorithmMultiCore::LocalCriticalValueList AlgorithmMultiCore::localCriticalValues;
AlgorithmMultiCore::LocalSPT AlgorithmMultiCore::spt;

AlgorithmMultiCore::AlgorithmMultiCore()
:   AlgorithmSingleCore(),
    c_diagonals_vector()
{
    auto copyspt = [&] () {
        return new PolygonShortestPathsFS(*qtri.local());
    };
    spt = LocalSPT(copyspt);
}


AlgorithmMultiCore::~AlgorithmMultiCore() 
{
    clear_qtri();
    localCriticalValues.clear();
    c_diagonals_vector.clear();
}

void AlgorithmMultiCore::clear_qtri()
{
    for (LocalTriangulation::iterator i = qtri.begin(); i != qtri.end(); ++i)
    {
        Triangulation* copy = *i;
        Q_ASSERT(copy != &Q->triangulation);    //  DO NOT delete the original
        delete copy;
    }
    for (LocalSPT::iterator i = spt.begin(); i != spt.end(); ++i)
    {
        PolygonShortestPathsFS* psf = *i;
        delete psf;
    }
    qtri.clear();
}

void AlgorithmMultiCore::triangulate()
{
    Algorithm::triangulate();

    //  local copies for Q->triangulation
    clear_qtri();
    qtri = LocalTriangulation([this] () -> frechet::poly::Triangulation* { 
        return new Triangulation(Q->triangulation); 
    });

    //  vector of c-diagonals
    c_diagonals_vector.clear();
    c_diagonals_vector.reserve(c_diagonals().size());
    std::copy(c_diagonals().begin(), c_diagonals().end(),
            std::back_inserter(c_diagonals_vector));
}


//
//  Calculate Reachability Graphs
//
void AlgorithmSingleCore::calculateReachabilityGraphs()
{
//    std::cout << "working in thread " << ConcurrencyContext::instance.currentThread() << std::endl;
    for(int i=0; i < l(); ++i) {
        testCancelFlag();   //  throws InterruptedException
        CRG(i, i+1);
    }
    Q_ASSERT(assertReachabilityGraphs());
}

void AlgorithmMultiCore::calculateReachabilityGraphs()
{
    auto lambda = [&](size_t i) {
        CRG(i, i+1);
        testCancelFlag();   //  throws InterruptedException
    };

    tbb::static_partitioner sp;
    tbb::parallel_for((size_t)0, (size_t)l(), lambda, sp);
    Q_ASSERT(assertReachabilityGraphs());
    //printHistogram();
}

bool AlgorithmSingleCore::assertReachabilityGraphs()
{
    for(int i=0; i < l(); ++i)
    {
        Triangulation::Edge e;
        GraphPtr rg1 = CRG(i-1,i,e);
        GraphPtr rg2 = CRG(i,i+1,e);
        GraphPtr rg3 = CRG(i+1,i+2,e);
        Q_ASSERT(Graph::is_adjacent_to(*rg1,*rg2));
        Q_ASSERT(Graph::is_adjacent_to(*rg2,*rg3));
    }
    return true;
}


void AlgorithmSingleCore::calculateValidPlacements(double epsilon)
{
    //  validPlacements contains only VERTICAL-VERTICAL edges.
    //  No horizontal edges.
    _validPlacements.clear();

    //  (1) finde die zu di und dj gehörigen freien Intervalle (aus dem Free-Space)
    //  (2) wähle aus jedem Intervall einen beliebigen Punkt: w1,...,wm
    //  wenn möglich, wählen wir Eckpunkte aus, sonst den kleinsten Punkt im Intervall
    _placements.resize(l());
    for(int i=0; i < l(); ++i)
        calculatePlacements(d(i), _placements[i]);

    Triangulation& qtri = Q->triangulation;
    PolygonShortestPathsFS spt(qtri);

    /*  Steps 7-11 */
    for(const Segment& d : c_diagonals())   //  parallel
    {
        //  for all placements in the free space
        int di = d.first;
        int dj = d.second;

        //  TODO calculate in one go, using mirrored R intervals ?!
        Algorithm::calculateValidPlacements(spt, epsilon, di, dj);
        Algorithm::calculateValidPlacements(spt, epsilon, dj, di);
    }
}

void AlgorithmMultiCore::calculateValidPlacements(double epsilon)
{
    //  validPlacements contains only VERTICAL-VERTICAL edges.
    //  No horizontal edges.
    _validPlacements.clear();

    //  (1) finde die zu di und dj gehörigen freien Intervalle (aus dem Free-Space)
    //  (2) wähle aus jedem Intervall einen beliebigen Punkt: w1,...,wm
    //  wenn möglich, wählen wir Eckpunkte aus, sonst den kleinsten Punkt im Intervall
    _placements.resize(l());

    //for(int i=0; i < l(); ++i)
    tbb::affinity_partitioner ap;
    auto lambda1 = [&] (size_t i) {
        calculatePlacements(d(i), _placements[i]);
    };
    tbb::parallel_for(0,l(),lambda1, ap);

    /*  Steps 7-11 */

    auto lambda2 = [&] (size_t i) {
        const Segment& d = c_diagonals_vector[i];
        int di = d.first;
        int dj = d.second;

        PolygonShortestPathsFS* local_spt = spt.local();

        //  in multi-threading mode, each thread holds a separate copy of Q->triangulation
/*      std::cout << "calculateValidPlacements "
                    << " thread " << ConcurrencyContext::instance.currentThread()
                    << " spt = " << &local_spt
                    << " tri = " << std::hex << &local_spt.triangulation()
                    << std::endl;
*/
#ifdef QT_DEBUG
        Triangulation& local_qtri = local_spt->triangulation();
        Q_ASSERT(&local_qtri != &Q->triangulation);
        Q_ASSERT(&local_qtri == this->qtri.local());
        Q_ASSERT(local_qtri.assertEquals(Q->triangulation));
#endif

        Algorithm::calculateValidPlacements(*local_spt, epsilon, di, dj);
        Algorithm::calculateValidPlacements(*local_spt, epsilon, dj, di);
    };

    tbb::parallel_for( (size_t)0, (size_t)c_diagonals_vector.size(), lambda2, ap);
}


//
//  The Main Loop
//

Graph::ptr AlgorithmSingleCore::CRG(int i, int j, Triangulation::Edge e)
{
    i = (i+l()) % l();
    j = (j+l()) % l();

    int h = mapKey(i,j);
    auto p = CRG_map.find(h);
    if (p!=CRG_map.end()) {
        return p->second;
    }
    else {
        Graph::ptr g = create_CRG(i,j,e);
        CRG_map.insert(std::make_pair(h,g));
        return g;
    }
}

Graph::ptr AlgorithmSingleCore::create_CRG(int i, int j, Triangulation::Edge e)
{
    Graph::ptr result;

#ifdef QT_DEBUG
    //  i,j are supposed to be incident to e
    if (e.first != Triangulation::Face_handle()) {
        Triangulation::VertexPair v2 = P->triangulation.incidentVertexes(e);
        if (v2.first->dindex == j)
            std::swap(v2.first, v2.second);
        Q_ASSERT(v2.first->dindex == i);
        Q_ASSERT(v2.second->dindex == j);
        Q_ASSERT(v2.first->pindex == d(i));
        Q_ASSERT(v2.second->pindex == d(j));
    }
    else if (e!=Triangulation::Edge()) {
        //  degenerated triangulation
        Q_ASSERT(l()==2);
        Q_ASSERT(i==0 && j==1 || i==1 && j==0);
    }
#endif

    if ((i+1)%l() == j%l()) {
        //  RG(i,i+1)
        Q_ASSERT(l()==2 || (e==Triangulation::Edge()) || Triangulation::is_polygon_edge(e));
        result = create_RG(i);
    }
    else {
        Q_ASSERT(Triangulation::is_inner_face(e.first));

        //      where (d_i,d_h,d_j) is a triangle in T_C'
        Triangulation::EdgePair e2 = P->triangulation.adjacentEdges(e,d(i),d(j));

        //std::cout << e << " " << e2.first << " " << e2.second << std::endl;

        Q_ASSERT( e2.first != e );
        Q_ASSERT( e2.second != e );
        //Q_ASSERT( P->triangulation.is_ccw_face(e2.first.first));
        //Q_ASSERT( P->triangulation.is_ccw_face(e2.second.first));
        Q_ASSERT(! P->triangulation.is_outer_edge(e2.first));
        Q_ASSERT(! P->triangulation.is_outer_edge(e2.second));

        Triangulation::Vertex_handle vh = e.first->vertex(e.second);
        Q_ASSERT(!Triangulation::is_south_pole(vh));

        Q_ASSERT(e2.first.first==e2.second.first);  //  same Face
        Q_ASSERT(e2.first.first->has_vertex(vh));
        Q_ASSERT(e2.second.first->has_vertex(vh));

        int h = vh->dindex;
        Q_ASSERT(h >= 0);
        Q_ASSERT(d(h)==vh->pindex);

        //  pick h so to avoid infinite recursion !

        P->triangulation.mirror(e2.first);
        Q_ASSERT(! P->triangulation.is_outer_edge(e2.first));

        //std::cout << e.first << std::endl;
        //std::cout << e2.second.first->neighbor(e2.second.second) << std::endl;
        P->triangulation.mirror(e2.second);
        Q_ASSERT(! P->triangulation.is_outer_edge(e2.second));

        //std::cout << "MERGE( CRG("<<i<<","<<h<<"), CRG("<<h<<","<<j<<"))" << std::endl;
        result = MERGE_inner( CRG(i,h,e2.first), CRG(h,j,e2.second) );
    }

    int di = d(i);
    int dj = d(j);
    if (is_c_diagonal(di,dj))
        COMBINE(result, di,dj);

    return result;
}


Graph::ptr AlgorithmSingleCore::create_RG(int i)
{
    return Algorithm::calculate_RG(i);
}


Graph::ptr AlgorithmSingleCore::MERGE_inner(const GraphPtr A, const GraphPtr B)
{
    //std::cout << "MERGE("<<A->origin.i<<","<<B->origin.i<<")" << std::endl;
    testCancelFlag();

    GraphPtr C = A->merge2(B.get());
    C->setOriginMerge2(A,B);
    return C;
}


Graph::ptr AlgorithmSingleCore::MERGE_outer(const GraphPtr A, const GraphPtr B)
{
    //  return A_HV * B_VV * A_VH
    GraphPtr C = A->merge3(B.get());
    C->setOriginMerge3(A,B);
    return C;

}


void AlgorithmSingleCore::COMBINE(GraphPtr A, int di, int dj)
{
    //  apply only to the VERTICAL-VERTICAL part of the matrix.
    GraphPtr P = getValidPlacements(di,dj);
    A->combine(P.get());
    A->setOriginCombine(P);
}


Graph::ptr AlgorithmSingleCore::findFeasiblePath()
{
    //FrechetViewApplication::instance()->pushTimer();

    Q_ASSERT(l() >= 2);
    for(int i=0; i < l(); ++i) //  parallel
    {
        testCancelFlag();

        Triangulation::Edge e = diagonalEdge(i);

        //  G = MERGE( RG(i-1,i), CRG(i,i-1,true), RG(i-1,i) );
        GraphPtr rg = CRG(i-1, i);
        GraphPtr G = MERGE_outer(rg, CRG(i, i-1, e));
        //  returns only the HH part of G

        //  Query G for a feasible path starting in [d_(i-1),d_i)
        //  (and ending in [d_(i-1) +n, d_i +n))
        Q_ASSERT(G->hmask() == rg->hmask());

        if (isSolution(G))
            return G;
    }

    //FrechetViewApplication::instance()->popTimer("Main Loop");

    return Graph::ptr();
}

bool AlgorithmSingleCore::isSolution(GraphPtr G)
{
    //  Note: search/foundDiagonalElement is decouopled for other implementations
    //  the default implementation returns a result immediately
    G->queryDiagonalElement();
    return (G->foundDiagonalElement() >= 0);
}


AlgorithmTopoSort::AlgorithmTopoSort() :
        AlgorithmMultiCore(),
        topo(), roots(), pool()
{ }

AlgorithmTopoSort::~AlgorithmTopoSort()
{ }

reach::Graph::ptr AlgorithmTopoSort::create_RG(int i) {
    //  calculate, enqueue - and insert into topo list
    Graph::ptr RG = Algorithm::calculate_RG(i);
    //  do not COMBINE
    Q_ASSERT(RG->origin.blevel == 0);
    addTopoSort(RG);
    return RG;
}

reach::Graph::ptr AlgorithmTopoSort::MERGE_inner(const reach::Graph::ptr A, const reach::Graph::ptr B) {
    //  DON'T calculate - just insert into topo list
    Graph::ptr C = newGraph(graphModel, graphModel->mergedHMask(A->hmask(),B->hmask()));
    C->setOriginMerge2(A, B);
    addTopoSort(C);
    return C;
}

reach::Graph::ptr AlgorithmTopoSort::MERGE_outer(const reach::Graph::ptr A, const reach::Graph::ptr B) {
    //  DON'T calculate - just insert into topo list
    GraphPtr C = newGraph(A->graphModel(), A->hmask());
    C->setOriginMerge3(A, B);
    addTopoSort(C);
    roots.push_back(C);
    return C;
}

void AlgorithmTopoSort::COMBINE(reach::Graph::ptr A, int di, int dj) {
    //  DON'T calculate - just update origin
    GraphPtr P = getValidPlacements(di, dj);
    A->setOriginCombine(P);
    P->origin.succesors++;
}

bool AlgorithmTopoSort::isSolution(GraphPtr G) {
    //  Solutions are evaluated later
    return false;
}

GraphPtr AlgorithmTopoSort::findFeasiblePath() {
    //  don't calculate - just insert Graphs into topo list
    AlgorithmMultiCore::findFeasiblePath();

    //  next: calculate from bottom upwards
    for (int blevel=0; blevel < topo.size(); ++blevel)
    {
        //  calculate one level
        for (Graph::ptr G : topo[blevel]) {
            calculate(G);
        }

        barrier(false); // before re-using matrix buffers

        for (Graph::ptr G : topo[blevel]) {
            G->resetConditions();
            // after a barrier, all conditions are met (we need not track them)
            reclaimPredecessors(G);
        }
    }
    //  collect resuts from roots
    barrier(true);
    for(Graph::ptr G : roots)
        if (G->foundDiagonalElement() >= 0)
            return G;
    return Graph::ptr();
}

void AlgorithmTopoSort::calculate(Graph::ptr G)
{
    Graph* A = G->origin.A.get();
    Graph* B = G->origin.B.get();
    Graph* P = G->origin.P.get();

    switch (G->origin.operation) {
        case Graph::Origin::RG:
            Q_ASSERT(!A && !B);
            break;
        case Graph::Origin::MERGE2:
            Q_ASSERT(A && B);
            G->merge2(A,B, &pool);
            break;

        case Graph::Origin::MERGE3:
            Q_ASSERT(A && B);
            Q_ASSERT(!P);
            G->merge3(A,B, &pool);
            G->queryDiagonalElement();
            break;
    }

    if (P)
        G->combine(P);
}

void AlgorithmTopoSort::reclaimPredecessors(Graph::ptr G)
{
    Graph* A = G->origin.A.get();
    Graph* B = G->origin.B.get();
    Graph* P = G->origin.P.get();

    if (A && (--A->origin.succesors==0)) // count dependencies on A_VV
        reclaim(A);
    if (B && (--B->origin.succesors==0)) // count dependencies on B_VV
        reclaim(B);
    if (P && (--P->origin.succesors==0))
        reclaim(P);
}

void AlgorithmTopoSort::reclaim(Graph* G)
{
    //  reclaim matrix data (will be used for upper levels)
    G->release(VERTICAL,VERTICAL,&pool);
    G->release(HORIZONTAL,HORIZONTAL,&pool);
    //  Note: G_HH may be used for finding a solution
    //  but the solution will be found before the next barrier()    
}

void AlgorithmTopoSort::cleanup()
{
    AlgorithmMultiCore::cleanup();
    topo.clear();
    roots.clear();
    pool.clear();
}

void AlgorithmTopoSort::addTopoSort(Graph::ptr G) {
    if (topo.size() <= G->origin.blevel)
        topo.resize(G->origin.blevel+1);
    topo[G->origin.blevel].push_back(G);
}

void AlgorithmTopoSort::barrier(bool finish) {
    if (ConcurrencyContext::hasGpuSupport())
    {
        if (finish)
            clFinish(ConcurrencyContext::clQueue());
        else
            clEnqueueBarrier(ConcurrencyContext::clQueue());
    }
}