kalgorithm.cpp

Go to the documentation of this file.

#include <k_frechet/kalgorithm.h>

using namespace frechet;
using namespace data;
using namespace k;
using namespace app;
using namespace fs;

kAlgorithm::kAlgorithm(FreeSpace::ptr fs)
    : components(fs->components()),
      domain(),
      working(), greedy(), bf()
{
    domain.H = Interval(0,fs->n-1);
    domain.V = Interval(0,fs->m-1);
    //setupMapsAndResults();
}


void kAlgorithm::setupMapsAndResults()
{
    Q_ASSERT(domain.H.lower()==0.0 && domain.V.lower()==0.0);

    BitSet::iterator i;
    size_t count = *components.last()+1;
    working.set = BitSet((int)count);
    greedy.result = BitSet((int)count);
    bf.result = BitSet((int)count);


    count = components.count();
    working.hmap.clear();
    working.vmap.clear();
    working.hmap.reserve(count);
    working.vmap.reserve(count);

    for(i = components.begin(); i.valid(); ++i)
    {
        //  i.offset() is the index in the Array = ID of component
        component_id_t compID = *i;
        const IntervalPair& ival = components.intervals() [compID];

        working.hmap.push_back(MappedInterval(ival.H,compID));
        working.vmap.push_back(MappedInterval(ival.V,compID));
    }

    std::sort(working.hmap.begin(),working.hmap.end());
    std::sort(working.vmap.begin(),working.vmap.end());

    greedy.k_x = greedy.k_y = greedy.k_xy = greedy.k_yx = UNKNOWN;
    bf.k_min = bf.k_max = bf.k_optimal = bf.k_current = bf.k_iteration = UNKNOWN;
}

int kAlgorithm::runGreedy()
{
    QMutexLocker guard(&working.lock);

    setupMapsAndResults();

    if (greedy_xy()==NO_SOLUTION) {
        greedy.k_yx =
        bf.k_min = bf.k_max = bf.k_optimal = bf.k_current = bf.k_iteration = NO_SOLUTION;

        emit greedyFinished();
        return NO_SOLUTION;
    }

    greedy.result.swap(working.set);
    greedy_yx();
    Q_ASSERT(greedy.k_xy>0 && greedy.k_yx>0);

    if (greedy.k_yx < greedy.k_xy)
        greedy.result.swap(working.set);

    /*  The results of the Greedy algorithm
     *  serve as lower and upper bounds for the Brute Force algorithm
     * */
    bf.k_min = greedy.lowerBound();
    bf.k_max = greedy.upperBound();

    if (bf.k_min == bf.k_max) {
        //  the Greedy result is already optimal. Fine.
        bf.k_optimal = bf.k_max;
        bf.result = greedy.result;
    }
    else {
        bf.k_optimal = UNKNOWN;
    }

    emit greedyFinished();
    return bf.k_max;
}


int kAlgorithm::greedy_xy()
{
    working.set.clear();

    greedy.k_x = greedy.k_xy = RUNNING;
    greedy.k_xy = greedy_12(
                working.hmap, domain.H,
                working.vmap, domain.V,
                greedy.k_x);
    Q_ASSERT(greedy.k_xy==NO_SOLUTION || greedy.k_xy==working.set.count());
    return greedy.k_xy;
}

int kAlgorithm::greedy_yx()
{
    working.set.clear();

    greedy.k_y = greedy.k_yx = RUNNING;
    greedy.k_yx = greedy_12(
                working.vmap, domain.V,
                working.hmap, domain.H,
                greedy.k_y);
    Q_ASSERT(greedy.k_yx==NO_SOLUTION || greedy.k_yx==working.set.count());
    return greedy.k_yx;
}

int kAlgorithm::greedy_12(
        const IntervalMap& map1, const Interval& range1,
        const IntervalMap& map2, const Interval& range2,
        int& k_1)
{
    k_1 = greedy_1(map1,range1);   //  first axis
    if (k_1 <= 0) return NO_SOLUTION;

    int k_2 = greedy_2(map2,range2);   //  second axis
    //  note: k_2 can be 0, if all intervals are already contained in k_1
    if (k_2 >= 0)
        return k_1 + k_2;
    else
        return NO_SOLUTION;
}


int kAlgorithm::greedy_1(const IntervalMap& map, const Interval& range)
{
    //  scan first axis
    int k = 0;
    int i = 0;
    double x = range.lower();

    while((x < range.upper()) && (i < map.size()))
    {   //  find the largest interval containing x
        if (map[i].lower() > x)
            return NO_SOLUTION;

        double x_next = x;
        component_id_t comp_next;

        for( ; (i < map.size()) && (map[i].lower() <= x); i++)
        {
            double xupper = map[i].upper();
            if (xupper > x_next) {
                x_next = xupper;
                comp_next = map[i].componentID;
            }
        }

        if (x_next > x) {
            x = x_next;
            working_add(comp_next);
            ++k;
        }
        else {
           return NO_SOLUTION;   //  we have made no progress
        }
    }

    if (x < range.upper())
        return NO_SOLUTION;
    else
        return k;
}

int kAlgorithm::greedy_2(const IntervalMap& map, const Interval& range)
{
    //  Scan second axis
    //  in addition, we must look for already selected intervals
    int k = 0;
    int i = 0;
    double y = range.lower();

    while((y < range.upper()) && (i < map.size()))
    {
        if (map[i].lower() > y)
            return NO_SOLUTION; //  a gap

        double y0 = y;
        double y_next = y;
        size_t comp_next;

        for( ; (i < map.size()) && (map[i].lower() <= y); ++i)
        {
            double yupper = map[i].upper();
            if (working.set.contains(map[i].componentID))
            {
                if (yupper > y)
                    y = yupper;
            }
            else if (yupper > y_next) {
                y_next = yupper;
                comp_next = map[i].componentID;
            }
        }

        if (y_next > y) {
            y = y_next;
            working_add(comp_next);
            ++k;
        }
        else if (y == y0) {
            return NO_SOLUTION;   //  we have made no progress
        }
    }

    if (y < range.upper())
        return NO_SOLUTION;
    else
        return k;
}

int kAlgorithm::runBruteForce(volatile bool* cancelFlag) throw (InterruptedException)
{
    QMutexLocker guard(&working.lock);

    if (bf.k_max == NO_SOLUTION)
        return NO_SOLUTION; //  Greedy was unsucessful. No reason to continue.

    Q_ASSERT(bf.k_min > 0 && bf.k_max > 0);
    //  these must have been calculated by the Greedy algorithm.

    if (bf.k_min == bf.k_max) {
        //  Greedy is already optimal. That was easy.
        bf.k_current = bf.k_iteration = bf.k_optimal = bf.k_max;
        bf.result = greedy.result;
        return bf.k_optimal;
    }

    //  TODO delegate to thread and return RUNNING immediately.
    bf.stack1.reserve(bf.k_max);
    bf.stack2.reserve(bf.k_max);

    bf.k_optimal = RUNNING;
    for(bf.k_iteration = bf.k_min; bf.k_iteration <= bf.k_max; bf.k_iteration++)
    {
        //  TODO parallel: we could run two searches in parallel (x-y and y-x)
        //  the fastest one wins
        //  for start, search only one way. (or: change direction every iteration?)
        //working.clear();
        //bf.k_current=0;

        emit bruteForceIteration(bf.k_iteration);
        if (bf_search_1(cancelFlag)) {
            //  found a result
            Q_ASSERT((bf.stack1.size()+bf.stack2.size()) == bf.k_iteration);
            //  copy stack2 to working set
            for( ; !bf.stack2.empty(); bf.stack2.pop_back()) {
                component_id_t comp = working.vmap[bf_top2().first].componentID;
                Q_ASSERT(!working.set.contains(comp));
                working.set += comp;
            }
            //  move working set to bf.result
            bf.result.swap(working.set);
            bf.k_optimal = bf.k_iteration;
            Q_ASSERT(bf.result.count() == bf.k_optimal);

            emit bruteForceFinished();
            return bf.k_optimal;
        }
    }
    //  We MUST have a result !!
    Q_ASSERT(false);
    return NO_SOLUTION;
}

void kAlgorithm::working_add(const MappedInterval &ival) {
    working_add(ival.componentID);
}

void kAlgorithm::working_remove(const MappedInterval &ival) {
    working_remove(ival.componentID);
}

void kAlgorithm::working_add(component_id_t comp)
{
    Q_ASSERT(!working.set.contains(comp));
    working.set += comp;
}

void kAlgorithm::working_remove(component_id_t comp)
{
    Q_ASSERT(working.set.contains(comp));
    working.set -= comp;
}

std::pair<int, double> &kAlgorithm::bf_top1() {
    Q_ASSERT(!bf.stack1.empty());
    return bf.stack1[bf.stack1.size()-1];
}

std::pair<int, double> &kAlgorithm::bf_top2() {
    Q_ASSERT(!bf.stack1.empty());
    return bf.stack2[bf.stack2.size()-1];
}

void kAlgorithm::bf_push1(int i, double x) {
    bf.stack1.push_back(std::make_pair(i,x));
    Q_ASSERT(working.set.count() == (bf.stack1.size()-1));
}

bool kAlgorithm::bf_pop1() {
    Q_ASSERT(bf.stack1.size() > 0);
    bf.stack1.pop_back();
    return !bf.stack1.empty();
}

void kAlgorithm::bf_push2(int i, double x) {
    bf.stack2.push_back(std::make_pair(i,x));
    //  don't maintain the working set
}

bool kAlgorithm::bf_backtrack2() {
    bf.stack2.clear();
    return false;
}

bool kAlgorithm::bf_search_1(volatile bool* cancelFlag) throw (InterruptedException)
{
    //  search horizontal range first
    const IntervalMap& map = working.hmap;
    Interval range = domain.H;

    bf.stack1.clear();
    bf.stack2.clear();
    working.set.clear();
    bf.stack1.push_back(std::make_pair(0,domain.H.lower()));

    for(;;)
    {
outer_loop:
        if (cancelFlag && *cancelFlag)
        {
            bf.k_current = bf.k_iteration = bf.k_optimal = UNKNOWN;
            throw InterruptedException();
        }

        Q_ASSERT(!bf.stack1.empty());
        int& i = bf_top1().first;
        double& x = bf_top1().second;

        for( ; (i < map.size()) && (map[i].lower() <= x); ++i)
        {
            if (map[i].upper() <= x) continue;  //  already covered, not interesting
            //  TODO sort choices?
            //  By interval length -> the first solution is the greedy solution
            //  By partial solution (stored from previous iteration) ?

            //  recurse
            if (map[i].upper() >= range.upper()) {
                working_add(map[i]);

                if (bf_search_2()) return true; //  switch to vertical

                working_remove(map[i]);
            }
            else if (bf.stack1.size() < bf.k_iteration) {
                //  push
                working_add(map[i]);
                bf_push1(i+1,map[i].upper());
                goto outer_loop;
            }
        }

        //  pop
        if (!bf_pop1()) return false;
        working_remove(map[bf_top1().first]);
        //  continue
        bf_top1().first++;
    }
    //  never come here
    Q_ASSERT(false);
    return false;
}

bool kAlgorithm::bf_search_2() throw (InterruptedException)
{    
//    if (abortRequested)
//        throw InterruptedException();

    const IntervalMap& map = working.vmap;
    Interval range = domain.V;

    //  Scan second axis
    int i = 0;
    double y = range.lower();

    while((y < range.upper()) && (i < map.size()))
    {
        if (map[i].lower() > y)
            return bf_backtrack2(); //  a gap

        double y0 = y;
        double y_next = y;
        int i_next = -1;

        for( ; (i < map.size()) && (map[i].lower() <= y); ++i)
        {
            double yupper = map[i].upper();
            if (working.set.contains(map[i].componentID))
            {
                if (yupper > y)
                    y = yupper;
            }
            else if (yupper > y_next) {
                y_next = yupper;
                i_next = i;
            }
        }

        if (y_next > y) {
            if ((bf.stack1.size()+bf.stack2.size()) >= bf.k_iteration)
                return bf_backtrack2(); //  reached k-limit

            Q_ASSERT(i_next >= 0);
            Q_ASSERT(!working.set.contains(map[i_next].componentID));
            bf_push2(i_next,y_next);
            y = y_next;
        }
        else if (y == y0)
            return bf_backtrack2();  //  we have made no progress
    }

    if (y < range.upper())
        return bf_backtrack2();
    else
        return true;
}


#define SORT_TIDY

bool MappedInterval::operator<(const MappedInterval &that) const {
    if (lower() != that.lower())
        return lower() < that.lower();
#ifdef SORT_TIDY
    if (upper() != that.upper())
        return upper() > that.upper();
    return componentID < that.componentID;
#else
    //  sort arbitrarily (but deterministic, of course)
    std::hash<const MappedInterval*> hsh;
    return hsh(this) < hsh(&that);
#endif
}

bool MappedInterval::operator==(const MappedInterval &that) const {
    return (lower()==that.lower())
            && (upper()==that.upper())
            && (componentID==that.componentID);
}