optimise.cpp

Go to the documentation of this file.

#include <poly/algorithm.h>
#include <poly/parallel.h>

#include <poly/poly_utils.h>
#include <frechet/poly/algorithm.h>
#include <data/numeric.h>
#include <app/concurrency.h>

#include <tbb/tbb.h>

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

#include <optimise_impl.h>

void Algorithm::collectCriticalValue(const accurate_float& x)
{
    if (x >= 0.0) {
        //  Aufrunden an der allerletzen Stelle.
        double dx = round_up((double)x,16);
        currentCriticalValueList().push_back(dx);
    }
}

Algorithm::CriticalValueList& Algorithm::currentCriticalValueList() {
    return criticalValues;
}

Algorithm::CriticalValueList& AlgorithmMultiCore::currentCriticalValueList() {
    return localCriticalValues.local();
}

void Algorithm::collectCriticalValues(bool a, bool b, bool c, bool d)
{
    const Curve& P = this->P->curve;
    const Curve& Q = this->Q->curve;

    //  (1) calculate critical values
    criticalValues.clear();
    if (a) collectCriticalValues_a();
    if (b) {
        collectCriticalValues_b(P,Q);
        collectCriticalValues_b(Q,P);
    }
    if (c) {
        collectCriticalValues_c(P,Q);
        collectCriticalValues_c(Q,P);
    }
    if (d) collectCriticalValues_d();

    testCancelFlag();
    combineCriticalValues();
    removeDuplicates(criticalValues);

    testCancelFlag();
    std::cout << "    Data:         "<<criticalValues.size()<<" critical values." << std::endl;
    printDebugTimer("Critical Values");
}


void AlgorithmSingleCore::combineCriticalValues()
{
    //  sort
    std::sort(criticalValues.begin(), criticalValues.end());
}

void AlgorithmMultiCore::combineCriticalValues()
{
    //  reduce localValueSets to one globalValueSet
    criticalValues = localCriticalValues.combine(AlgorithmMultiCore::combine);
    //  sort
    tbb::parallel_sort(criticalValues.begin(),criticalValues.end());

    localCriticalValues.clear();
}


void Algorithm::removeDuplicates(CriticalValueList& vector)
{
    double* begin = &vector[0];
    double* end = begin+vector.size();
    double* a = begin;
    double* b = begin;

    double eps = round_up(0.0,16);
    while(a<end) {
        while(*a-*b < eps)
            if (++a>=end) break;
        *++b = *a++;
    }

    vector.resize(b+1 - begin);
}


const AlgorithmMultiCore::CriticalValueList& AlgorithmMultiCore::combine(
        CriticalValueList& a, const CriticalValueList& b)
{
    size_t s1 = a.size();
    size_t s2 = b.size();
    a.resize(s1+s2);
    memcpy(&a[s1], &b[0], s2*sizeof(double));
    //b.clear();
    return a;
}

void Algorithm::collectCriticalValues_a()
{
    const Curve& P = this->P->curve;
    const Curve& Q = this->Q->curve;

    //  Critical Values Type (a)
    //  distance between endpoints; applies only to open curves
    if (!P.isClosed() && !Q.isClosed())
    {
        collectCriticalValue( euclidean_distance<accurate_float>(P.front(),Q.front()));
        collectCriticalValue( euclidean_distance<accurate_float>(P.back(),Q.back()));
    }
}

void Algorithm::collectCriticalValues_b(const Point& p, const Curve& Q)
{
    //std::cout << " Thread "<<ConcurrencyContext::instance.currentThread()<<" collect b ["<<i<<".."<<j<<"]"<<std::endl;
    for(int k=1; k < Q.size(); ++k)
    {
        QLineF segment (Q[k-1],Q[k]);
        accurate_float dist = segment_distance<accurate_float>(segment,p);
        collectCriticalValue(dist);
    }
    testCancelFlag();
}

void AlgorithmSingleCore::collectCriticalValues_b(const Curve& P, const Curve& Q)
{
    //  Critical Values Type (b)
    //  distance from a point to a segment
    for(int i=0; i < P.length(); ++i)
        Algorithm::collectCriticalValues_b(P[i], Q);
}

void AlgorithmMultiCore::collectCriticalValues_b(const Curve& P, const Curve& Q)
{
    //  Critical Values Type (b)
    //  distance from a point to a segment
    auto lambda = [&](size_t i) {
        Algorithm::collectCriticalValues_b(P[i], Q);
    };
    //  TODO [] (size_t i)

    tbb::parallel_for( (size_t)0, (size_t)P.size(), lambda);
}

void Algorithm::collectCriticalValues_c(const QLineF& segment, const Curve& Q)
{
    //  for each segment of P
    //      for each pair of vertices in Q
    //          distance to intersection between segment
    //          and bisector
    //std::cout << " Thread "<<ConcurrencyContext::instance.currentThread()<<" collect c ["<<i<<".."<<j<<"]"<<std::endl;

    for(int k=0; k < Q.size(); ++k) {
        for (int l = k+1; l < Q.size(); ++l) {
            /*  detecting intersections: allow some tolerance.
             *  Better find too many critical values than missing one !!
             */
            accurate_float dist = bisector_distance<accurate_float> (Q[k], Q[l], segment, 1e-3);
            collectCriticalValue(dist);
        }
        testCancelFlag();
    }
}

void AlgorithmSingleCore::collectCriticalValues_c(const Curve& P, const Curve& Q)
{
    for(int i=1; i < P.length(); ++i) {
        QLineF segment(P[i-1],P[i]);
        Algorithm::collectCriticalValues_c(segment, Q);
    }
}

void AlgorithmMultiCore::collectCriticalValues_c(const Curve& P, const Curve& Q)
{
    auto lambda = [&](size_t i) {
        QLineF segment(P[i-1],P[i]);
        Algorithm::collectCriticalValues_c(segment, Q);
    };

    tbb::parallel_for( (size_t)1, (size_t)P.length(), lambda);
}

void Algorithm::collectCriticalValues_d(const QLineF& diagonal)
{
    //  Critical values type (d)
    //  between a diagonal and reflex vertices
//    const Curve& P = this->P->curve;
    const Curve& Q = this->Q->curve;
    const Polygon& R = this->Q->reflex;
    Q_ASSERT(R.size() >= 1); //  da Q nicht konvex ist

    //std::cout << " Thread "<<ConcurrencyContext::instance.currentThread()<<" collect d ["<<r1<<"]"<<std::endl;

    for(int r1=0; r1 < R.size(); ++r1)
    {
        const Point& q1 = Q[R[r1]];
        //  type (b)
        accurate_float dist = segment_distance<accurate_float>(diagonal,q1);
        collectCriticalValue(dist);
        //  type (c)
        for (int r2 = r1+1; r2 < R.size(); ++r2) {
            const Point& q2 = Q[R[r2]];
            /*  detecting intersections: allow some tolerance.
             *  Better find too many critical values than missing one !!
             */
            dist = bisector_distance<accurate_float> (q1, q2, diagonal, 1e-3);
            collectCriticalValue(dist);
        }
        testCancelFlag();
    }
}

//for(const Segment& d : c_diagonals())
//QLineF segment(P[d.first], P[d.second]);

void AlgorithmSingleCore::collectCriticalValues_d()
{
    const Curve& P = this->P->curve;
    for(const Segment& d : c_diagonals())
    {
        QLineF diagonal(P[d.first], P[d.second]);
        Algorithm::collectCriticalValues_d(diagonal);
    }
}

void AlgorithmMultiCore::collectCriticalValues_d()
{
    const Curve& P = this->P->curve;
    auto lambda = [&] (size_t i) {
        Segment d = c_diagonals_vector[i];
        QLineF diagonal(P[d.first], P[d.second]);
        Algorithm::collectCriticalValues_d(diagonal);
    };

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



double Algorithm::optimisePolygon(frechet::reach::FSPath::ptr path, double approximation)
{
    Q_ASSERT(canOptimisePoly());

    cleanup();
    setStatus((approximation==0.0) ? RUNNING_OPT_POLY : RUNNING_APPROX_POLY);

    GraphPtr result;
    double epsilon;
    Pointer ignored;
    reach::FSPath::ptr no_path;

    local_fs.reset(new FreeSpace(P->curve, Q->curve));

    if (approximation==0.0)
    {        
        //  === Binary Search on Critical Values ===

        //  (1) calculate critical values
        collectCriticalValues(false,true,true,true);

        // (2) binsearch on critical values
        int lo = -1;
        int hi = criticalValues.size();

        //      find lower bound through decideCurve
//      find higher bound through leaping bin search
        lo = binarySearch<Pointer> (criticalValues, lo, hi,
                                    &Algorithm::decideCurve, ignored);
        Q_ASSERT(lo >= 0 && lo < criticalValues.size());
        printDebugTimer("Lower Bound");

        epsilon = criticalValues[lo];
        local_fs->calculateFreeSpace(epsilon);
        result = decidePolygon(local_fs, no_path, epsilon, false);

        if (result) goto label_result;
        //  find hi
        int step = 1;           
        for (;;) {
            hi = lo + step;
            if (hi >= criticalValues.size()) {
                hi = criticalValues.size();
                break;
            }

            epsilon = criticalValues[hi];
            local_fs->calculateFreeSpace(epsilon);
            result = decidePolygon(local_fs, no_path, epsilon, false);

            if (result && step == 1) goto label_result;
            if (result)
                break;  //  found hi
            else {
                lo = hi;
                step *= 2;
            }
        }

        int i = binarySearch(criticalValues, lo, hi, &Algorithm::decidePolygon, result);
        Q_ASSERT(i >= 0 && i < criticalValues.size());
        epsilon = criticalValues[i];
    }
    else
    {
        //  === Nested Interval search ===
        double upper_bound = frechet::max_euclidean_distance(P->curve, Q->curve) + 1.0;
        //  d_F for Curve constitues the lower bound (is quicker to calculate!)
        epsilon = intervalSearch(0.0, upper_bound, approximation, &Algorithm::decideCurve, ignored);
        local_fs->calculateFreeSpace(epsilon);
        result = decidePolygon(local_fs, no_path, epsilon, false);

        if (result) goto label_result;
        epsilon = intervalSearch<GraphPtr> (epsilon,upper_bound, approximation, &Algorithm::decidePolygon, result);
    }

label_result:
    Q_ASSERT((bool)result);
    local_fs.reset();

    if (path) {
        //  FSPath is very sensitive to round-off errors.
        //  Because it is only needed for visualisation, we can afford
        //  a bit of tolerance.
        //  TODO reduce round-off errors with float128.
        double repsilon = round_up(epsilon,32);
        updatePath(path, boost::static_pointer_cast<Graph>(result),repsilon);
    }

    setStatus( result ? YES:NO );

    if (path) {
        emit optimisationResult(epsilon);
        emit pathUpdated((bool)result); //  only after status has changed
    }

    this->fs.reset();
    criticalValues.clear();

    cleanup();
    return epsilon;
}

double Algorithm::optimiseCurve(double approximation)
{
    Q_ASSERT(canOptimiseCurve());

    setStatus((approximation==0.0) ? RUNNING_OPT_CURVE : RUNNING_APPROX_CURVE);

    reach::Pointer result;
    double epsilon;
    local_fs.reset(new FreeSpace(P->curve, Q->curve));

    if (approximation==0.0)
    {
        //  === Binary Search on Critical Values ===
        //  (1) calculate critical values
        collectCriticalValues(true,true,true,false);

        int i = binarySearch<reach::Pointer> (criticalValues, -1, criticalValues.size(),
                                                &Algorithm::decideCurve, result);
        Q_ASSERT(i >= 0 && i < criticalValues.size());
        epsilon = criticalValues[i];
    }
    else
    {
        //  === Nested Interval search ===
        double upper_bound = frechet::max_euclidean_distance(P->curve,Q->curve)+1.0;
        epsilon = intervalSearch<reach::Pointer> (0.0,upper_bound, approximation, &Algorithm::decideCurve, result);
    }

    Q_ASSERT((bool)result);
    local_fs.reset();
    criticalValues.clear();

    this->fs.reset();
    setStatus(result ? YES:NO);

    emit optimisationResult(epsilon);
    //return startPointer;
    return epsilon;
}