main page namespaces classes files
ELL home

src/ell/Basin.cc

Go to the documentation of this file.

#include "ell/Basin.hh"
#include <biu/assertbiu.hh>

//##############################################################################
//##############################################################################

namespace ell
{

    Basin::~Basin()
    {}


} // namespace ell

//##############################################################################
//##############################################################################

#include "ell/util/PriorityQueue.hh"
#include <list>

namespace ell {

    Basin_FME::QueueVal::QueueVal()
     :  minNeigh(), minNeighE(0.0), minNeighGWE(UINT_MAX)
    {
    }
    
    Basin_FME::QueueVal::~QueueVal()
    {
    }


    Basin_FME::Basin_FME(   const double maxE_, 
                    const double deltaE_, 
                    const size_t maxToStore_)
      : maxE(maxE_),
        deltaE(deltaE_),
        maxToStore(maxToStore_)
    {}
    
    Basin_FME::~Basin_FME()
    {}
    
    double
    Basin_FME::enumerate(   const State* const localMin, 
                            StateCollector* scBasin,
                            StatePairCollector* scSurface)
    {
        return Basin_FME::flood( localMin, scBasin, 
                                maxE, deltaE, maxToStore, scSurface);
    }

    double 
    Basin_FME::flood(   
                        const State* const localMin, 
                        StateCollector* scBasin,
                        const double maxE_, 
                        const double deltaE, 
                        const size_t maxToStore, 
                        StatePairCollector* scSurface
            ) 
    {
        assertbiu( localMin != NULL, "No local Min given! (==NULL)"); 
        assertbiu( scBasin != NULL, "No StateCollector given! (scBasin==NULL)"); 
        assertbiu( localMin->getEnergy() <= maxE_, "E(localMin) <= maxE");
        
        using namespace ell::util;
        
          // the priority queue to store neighbors to handle in
        PriorityQueue<QueueVal> pq;
          // init using the given local minimum
        PriorityQueue<QueueVal>::InsertResult it = pq.insert(localMin);
        it.first->second.minNeigh.clear(); // has no minimal neighbor (is min)
        it.first->second.minNeighE = it.first->first.E;
        it.first->second.minNeighGWE = 0;  
        
        double maxE = maxE_;
        
        typedef std::list< PriorityQueue<QueueVal>::QueueKey  > NeighList;
        
        State* topState = localMin->clone();
        
        while (!pq.empty()) {
              // get top element
            PriorityQueue<QueueVal>::const_iterator top = pq.top();
              // uncompress top element into new State
            topState->uncompress(top->first.s);
            CSequence minNeigh; // neighbor with minimal energy
            double minNeighE = (double)INT_MAX; // energy of the minimal neighbor
             
            
              // get neighbor list
            State::NeighborListPtr neighbors = topState->getNeighborList();
            State::NeighborList::Iterator it = neighbors->begin();
              // contains all neighbors with E(top) < E <= maxE
            NeighList toStore;
              // the compressed sequence of the current neighbor checked
            CSequence next;
              // iterate through all neighbors
            while (neighbors->end() != it) 
            {
                  // get neighbor data
                double nextE = it->getEnergy();
                next = it->compress( next );
                  // increase iterator for next iteration
                ++it;

                  // check for lowest neighbor
                if (    (nextE < minNeighE) ||
                        (nextE == minNeighE && next < minNeigh ) ) 
                {
                    minNeigh = next;
                    minNeighE = nextE;
                }
                  // neighbor rejection
                if (    nextE <= maxE   // energy does not exceed maxE
                        && ( nextE > top->first.E  // energy higher than top or
                        || (nextE == top->first.E && next > top->first.s) 
                                // greater in sequence order 
                        ) ) 
                {
                    // store neighbor with higher energy that has to be added to pq
                    toStore.push_front(NeighList::value_type(nextE, next));
                }
                    
            }
              // handle top depending on its basin membership           
            if (    (minNeighE == top->second.minNeighE 
                    && minNeigh == top->second.minNeigh)   // belongs to the basin
                    || top->second.minNeigh.empty() ) // top is the local minimum
            {
                  // store because top belongs to the basin 
                scBasin->add(*topState);
                
                PriorityQueue<QueueVal>::iterator toFill;
                  // sort list with increasing energy
                toStore.sort();
                for (   NeighList::const_iterator n = toStore.begin(); 
                        n != toStore.end() && n->E <= maxE; 
                        n++) 
                {
                    toFill = pq.find(*n);
                    if (toFill == pq.end()) {
                          // insert empty value object
                        toFill = pq.insert(*n).first;
                          // feed data
                        toFill->second.minNeigh = top->first.s;
                        toFill->second.minNeighE = top->first.E;
                        toFill->second.minNeighGWE = top->second.minNeighGWE;
                          // check queue size and reduce if necessary
                        if ( pq.size() >= maxToStore ) { // queue to big
                              // get current maximal E in queue
                            double curMaxE = pq.getMaxE();
                              // reduce maximally allowed energy below cur. max
                            do { maxE -= deltaE; } while (maxE > curMaxE);
                              // erase all elements > maxE from the queue;
                            pq.reduceMaxE(maxE+0.00000001);  // <-- this is a HACK
                        }
                    } else  // update of entry if this one is smaller
                        if (    toFill->second.minNeighE > top->first.E ||
                                ( toFill->second.minNeighE == top->first.E &&
                                  toFill->second.minNeigh > top->first.s) )
                    {
                          // feed data
                        toFill->second.minNeigh = top->first.s;
                        toFill->second.minNeighE = top->first.E;
                        toFill->second.minNeighGWE = top->second.minNeighGWE;
                    }
                }
                
//          } else { // belongs NOT to the basin but is surface
                
                // UNFORTUNATELY WE DONT KNOW WHAT STATES OF THE BASIN ARE 
                // NEIGHBORED TO THIS STATE --> SO NO REPORT TO THE SURFACE
                // COLLECTOR IS DONE
            } 

              // remove top element
            pq.pop();           
        }
          // clear memory
        delete topState;
         
        return maxE;
    }

} // namespace ell

//##############################################################################
//##############################################################################


#include "ell/util/PriorityQueue.hh"
#include <list>
#include <set>

namespace ell {

    Basin_F::QueueVal::QueueVal()
     :  minNeigh(), minNeighE(0.0), minNeighGWE(UINT_MAX)
    {
    }
    
    Basin_F::QueueVal::~QueueVal()
    {
    }


    Basin_F::Basin_F(   const double maxE_, 
                        const double deltaE_, 
                        const size_t maxToStore_ )
      : maxE(maxE_),
        deltaE(deltaE_),
        maxToStore(maxToStore_)
    {}
    
    Basin_F::~Basin_F()
    {}
    
    double
    Basin_F::enumerate( const State* const localMin, 
                            StateCollector* scBasin,
                            StatePairCollector* scSurface)
    {
        return Basin_F::flood( localMin, scBasin, 
                                maxE, deltaE, maxToStore, scSurface);
    }

    double 
    Basin_F::flood( 
                        const State* const localMin, 
                        StateCollector* scBasin,
                        const double maxE_, 
                        const double deltaE, 
                        const size_t maxToStore, 
                        StatePairCollector* scSurface
            ) 
    {
        assertbiu( localMin != NULL, "No local Min given! (==NULL)"); 
        assertbiu( scBasin != NULL, "No StateCollector given! (scBasin==NULL)"); 
        assertbiu( localMin->getEnergy() <= maxE_, "E(localMin) <= maxE");
        
        using namespace ell::util;
        
          // the priority queue to store neighbors to handle in
        PriorityQueue<QueueVal> pq;
          // the hash to store seen elements of the basin
        typedef std::set<CSequence> BasinHash;
        BasinHash handled;
          // init using the given local minimum
        PriorityQueue<QueueVal>::InsertResult it = pq.insert(localMin);
        it.first->second.minNeigh.clear(); // has no minimal neighbor (is min)
        it.first->second.minNeighE = it.first->first.E;
        it.first->second.minNeighGWE = 0;  
        
        double maxE = maxE_;
        
        typedef std::list< PriorityQueue<QueueVal>::QueueKey  > NeighList;
        
          // create temporary objects that are altered during the flooding
        State* topState = localMin->clone();
        State* curNeigh = localMin->clone();
        
        while (!pq.empty()) {
              // get top element
            PriorityQueue<QueueVal>::const_iterator top = pq.top();
              // uncompress top element into new State
            topState->uncompress(top->first.s);
            CSequence minNeigh; // neighbor with minimal energy
            double minNeighE = (double)INT_MAX; // energy of the minimal neighbor
             
            
              // get neighbor list
            State::NeighborListPtr neighbors = topState->getNeighborList();
            State::NeighborList::Iterator it = neighbors->begin();
              // contains all neighbors with E(top) < E <= maxE
            NeighList toStore;
              // contains all neighbors with E < E(top) 
            NeighList toCheck;
              // the compressed sequence of the current neighbor checked
            CSequence next;
              // iterate through all neighbors
            while (neighbors->end() != it) 
            {
                  // get neighbor data
                double nextE = it->getEnergy();
                next = it->compress( next );
                  // increase iterator for next iteration
                ++it;

                  // check for lowest neighbor
                if (    (nextE < minNeighE) ||
                        (nextE == minNeighE && next < minNeigh ) ) 
                {
                    minNeigh = next;
                    minNeighE = nextE;
                }
                  // check if neighbor is of interest and not rejected
                if ( nextE <= maxE )    // energy does not exceed maxE
                {
                      // check if neighbor is of higher order --> to PQ
                    if (    nextE > top->first.E  // energy higher than top or
                            || (nextE == top->first.E && next > top->first.s) ) 
                    {        // greater in sequence order
                        // store neighbor with higher energy that has to be added to pq
                        toStore.push_front(NeighList::value_type(nextE, next));
                        
                    } else if (scSurface != NULL) {
                        // store all neighbors with lower energy 
                        toCheck.push_back(NeighList::value_type(nextE, next));
                    }
                }
                
            }
              // handle top depending on its basin membership           
            if (    (minNeighE == top->second.minNeighE 
                    && minNeigh == top->second.minNeigh)   // belongs to the basin
                    || top->second.minNeigh.empty() ) // top is the local minimum
            {
                  // store because top belongs to the basin 
                scBasin->add(*topState);
                  // add to hashed elements and mark that part of the basin
                handled.insert(top->first.s);
                
                  // put all neighbors with higher energy to PQ
                PriorityQueue<QueueVal>::iterator toFill;
                  // sort list with increasing energy
                toStore.sort();
                for (   NeighList::const_iterator n = toStore.begin(); 
                        n != toStore.end() && n->E <= maxE; 
                        n++) 
                {
                    toFill = pq.find(*n);
                    if (toFill == pq.end()) {
                          // insert empty value object
                        toFill = pq.insert(*n).first;
                          // feed data
                        toFill->second.minNeigh = top->first.s;
                        toFill->second.minNeighE = top->first.E;
                        toFill->second.minNeighGWE = top->second.minNeighGWE;
                          // check queue size and reduce if necessary
                        if ( pq.size() >= maxToStore ) { // queue to big
                              // get current maximal E in queue
                            double curMaxE = pq.getMaxE();
                              // reduce maximally allowed energy below cur. max
                            do { maxE -= deltaE; } while (maxE > curMaxE);
                              // erase all elements >= maxE from the queue;
                            pq.reduceMaxE(maxE);
                        }
                    } else  // update of entry if this one is smaller
                        if (    toFill->second.minNeighE > top->first.E ||
                                ( toFill->second.minNeighE == top->first.E &&
                                  toFill->second.minNeigh > top->first.s) )
                    {
                          // feed data
                        toFill->second.minNeigh = top->first.s;
                        toFill->second.minNeighE = top->first.E;
                        toFill->second.minNeighGWE = top->second.minNeighGWE;
                    }
                }
                  // check if surface / contact plane is of interest
                if (scSurface != NULL) {
                      // check all neighbors with lower energy if NOT part of the
                      // basin to enumerate the surface of the basin
                    for (   NeighList::const_iterator n = toCheck.begin(); 
                            n != toCheck.end(); n++) 
                    {
                          // check if neighbor is NOT part of basin
                        if ( handled.find(n->s) == handled.end() ) {
                              // uncompress neighbor
                            curNeigh->uncompress(n->s);
                              // add to surface reporter
                            scSurface->add( topState, curNeigh );
                        }
                    }
                }
                
            } else { // belongs NOT to the basin but in surface
                
                  // check if surface / contact plane is of interest
                if (scSurface != NULL) {
                      // check all neighbors with lower energy IF part of the
                      // basin to enumerate the surface of the basin
                    for (   NeighList::const_iterator n = toCheck.begin(); 
                            n != toCheck.end(); n++) 
                    {
                          // check if neighbor IS part of basin
                        if ( handled.find(n->s) != handled.end() ) {
                              // uncompress neighbor
                            curNeigh->uncompress(n->s);
                              // add to surface reporter
                            scSurface->add( curNeigh, topState );
                        }
                    }
                }
            } 

              // remove top element
            pq.pop();           
        }
          // clear memory
        delete topState;
        delete curNeigh;
         
        return maxE;
    }

} // namespace ell

//##############################################################################
//##############################################################################