main page namespaces classes files
ELL home

src/ell/LT_BarrierTree.cc

Go to the documentation of this file.

#include "ell/LT_BarrierTree.hh"
#include "ell/LT_SaddleNet.hh"

#include <iomanip>
#include <fstream>
#include <sstream>
#include <cmath>
#include <algorithm>
#include <list>
#include <stack>
#include <set>

#include <cstring>
#include <ctime>
    
namespace ell {

    const std::string LT_BarrierTree::LT_ID = "BT";


    LT_BarrierTree::LT_BarrierTree( const double maxMergeBarrier_ )
     :  LT_MinimaSet()
        , maxMergeBarrier(maxMergeBarrier_)
        , root(NULL)
        , min2node()
        , numOfLeaves(0)
    { 
    }

    LT_BarrierTree::LT_BarrierTree( const LT_BarrierTree& toCopy )
    {
        this->operator =(toCopy);
    }
    
    
    LT_BarrierTree&
    LT_BarrierTree::
	operator = ( const LT_MinimaSet& toCopy )
    {
          // check if toCopy is not the same object than this
        if (this == &toCopy) {
            return *this;
        }
          // check if it is a barrier tree
        const LT_BarrierTree* bt = dynamic_cast<const LT_BarrierTree*>(&toCopy);
        
        if (bt != NULL) {
            return this->operator =(*bt);
        } else {
            double oldMaxMergeBarrier = maxMergeBarrier;
              // clear data structure
            this->clear();
              // reset maximal merge barrier to previous value
            maxMergeBarrier = oldMaxMergeBarrier;
              // copy all minima
            for (size_t i=0; i<toCopy.getMinCount(); i++) {
                  // add minimum to tree
                this->addMin(toCopy.getMin(i));
            }
              // update mfe index
            mfeIndex = getMinIndex(toCopy.getMFEState());
              // update sorted status
            if (!toCopy.isSorted()) {
                sorted = false;
            }
        }
        return *this;
    }
    
    LT_BarrierTree&
    LT_BarrierTree::
	operator = ( const LT_BarrierTree& toCopy )
    {
          // check if toCopy is not the same object than this
        if (this == &toCopy) {
            return *this;
        }
          // clear data structures
        this->clear();
        
          // copy maximal merge barrier used in the other tree
        maxMergeBarrier = toCopy.maxMergeBarrier;
        
          // check if the other tree contains nodes
        if (toCopy.root != NULL) {
            root = new Node(NULL,NULL,NULL,true);
        } else {
            return *this;
        }
        
          // copy all nodes and update vMinima list 
        typedef std::pair<const Node*,Node*> NP;
        std::list< NP > nodes;
        nodes.push_back(NP(toCopy.root,root));
        
        State* dummyState = toCopy.getMFEState()->clone();
        
        numOfLeaves = 0;
        while (!nodes.empty()) {
              // get next node pair to handle
            NP n = nodes.front();
            nodes.pop_front();
            
            if (n.first != NULL) {
                  // check if internal node
                if (n.first->children != NULL) {
                    // push all children nodes to node stack 
                    const NodeList* children = n.first->children;
                    n.second->children = new NodeList(root, NULL);
                    NodeList** myChildren = &(n.second->children);
                    while (children != NULL) {
                        assertbiu(children->current != NULL, "current child of toCopy node has no link to other node (==NULL)");
                          // add new child to this tree
                        Node* child = new Node(n.second,NULL,NULL,true);
                        (*myChildren)->current = child;
                          // push to stack for later filling
                        nodes.push_back(NP(children->current,child));
                          // go to next child in the list
                        children = children->next;
                        myChildren = &((*myChildren)->next);
                        if (children != NULL)
                            *myChildren = new NodeList(root, NULL);
                    }
                    assertbiu(n.second->children != NULL , "internal node has no children after copying");
                      // copy all barrier states
                    const StateList* states = n.first->content;
                    n.second->content = new StateList(dummyState,NULL);
                    StateList* myStates = n.second->content;
                    while (states != NULL) {
                        assertbiu(states->current != NULL, "current state of toCopy node has no link to other node (==NULL)");
                        myStates->current = states->current->clone();
                          // go to next barrier
                        states = states->next;
                        myStates = myStates->next;
                        if (states != NULL)
                            myStates = new StateList(dummyState,NULL);
                    }
                    assertbiu(n.second->content != NULL || n.second->parent == NULL, "internal node has no content after copying");
                } else {
                    // leaf node
                    numOfLeaves++;
                    
                      // copy all minimum states and add to vMinima list
                    const StateList* states = n.first->content;
                    n.second->content = new StateList(dummyState,NULL);
                    StateList* myStates = n.second->content;
                    while (states != NULL) {
                        assertbiu(states->current != NULL
                                , "current state of toCopy node has no link to other node (==NULL)");
                          // add to list
                        State* s = states->current->clone();
                        myStates->current = s;
                          // add to vMinima
                        vMinima.push_back(s);
                          // update mapping
                        min2node[s] = n.second;
                          // go to next minimum in the list
                        states = states->next;
                        if (states != NULL)
                            myStates->next = new StateList(dummyState,NULL);
                        myStates = myStates->next;
                    }
                    assertbiu(n.second->content != NULL 
                                , "leaf node has no content after copying");
                }
            }
        }  // while (on all nodes)
        delete dummyState;
        assertbiu(numOfLeaves == toCopy.numOfLeaves
                        , "different number of leaves after copying");
        
          // get mfe state index
        mfeIndex = getMinIndex( toCopy.getMFEState() );
        assertbiu(mfeIndex != INVALID_INDEX
                        , "mfe structure not found after copying");
        
        sorted = false;
          // ensure sorting if origin is sorted
        if (toCopy.sorted) {
            this->sort();
        }
        
        return *this;
    }


    void 
    LT_BarrierTree::clear() 
    
    {
        if (root != NULL && mfeIndex != LandscapeTopology::INVALID_INDEX ) {
            delete root;
        }
        vMinima.clear();
        mfeIndex = LandscapeTopology::INVALID_INDEX;
        sorted = true;
        maxMergeBarrier = 0.0;
        root = NULL;
        min2node.clear();
        numOfLeaves = 0;
    }

    LT_BarrierTree::~LT_BarrierTree() {
        clear();
    }

    bool
    LT_BarrierTree::
	addSaddle(const size_t i, const size_t j, const State* const s) 
    
    {
          // check
        assertbiu(i!=j, "index i and j are equal!");
        assertbiu(i<vMinima.size(), "index i out of vMinima range");
        assertbiu(j<vMinima.size(), "index j out of vMinima range");
        assertbiu(s->getEnergy() >= vMinima[i]->getEnergy()
                        , "saddle has lower energy than minimum i");
        assertbiu(s->getEnergy() >= vMinima[j]->getEnergy()
                        , "saddle has lower energy than minimum j");
        
          // get nodes of the minima
        Node* nm1 = min2node[ vMinima[i] ];
        Node* nm2 = min2node[ vMinima[j] ];
        
          // check if minima not already merged
        if (nm1 == nm2) {
            return false;
        }
        
          // get representive minimum of the nodes
        const State* m1 = nm1->content->current;
        const State* m2 = nm2->content->current;
        
          // m1 is the smaller 'leading' minimum of both 
        if (State::less(m2,m1)) {
            std::swap(m1,m2);
            std::swap(nm1,nm2);
        }
        
//{
//std::cerr <<"\n######  ADDING OF (" <<i <<"=" <<vMinima[i]->toString() <<", " <<j<<"=" <<vMinima[j]->toString()<<", "<<s->toString() <<") ###\n";
// bool nomore = false;
//for (std::map<const State*,Node*>::const_iterator it=min2node.begin();it!=min2node.end();it++) {
//  std::cerr <<"## " <<getMinIndex(it->first) <<" = " <<it->first->toString() <<" --> " <<it->second;
//  Node* p = it->second->parent;
//  size_t i =0;
//  while (p!=NULL) {
//      std::cerr <<" --> " <<p;
//      p = p->parent;
//      if (i++ > 20) {
//          nomore = true;
//          break;
//      }
//  } std::cerr <<"\n";
//}
//std::cerr <<"## nm1 = " <<nm1 <<" = " <<m1->toString() <<" p=" <<nm1->parent <<" -> " <<nm1->parent->parent<<"\n";
//std::cerr <<"## nm2 = " <<nm2 <<" = " <<m2->toString() <<" p=" <<nm2->parent <<" -> " <<nm2->parent->parent<<"\n";
//}
          // get energy difference of higher minimum to saddle
        const double eDiff = s->getEnergy() - m2->getEnergy();
          // check if m2 has to be merged with m1
        if ( eDiff <= maxMergeBarrier) 
        {

            // move list of minima that were already merged with m2
            {
                StateList* l = nm2->content;
                while (l != NULL) {
                      // update node index of all minima that are merged with m2
                    min2node[ l->current ] = nm1;
                      // add all minima that are merged with m2 to m1
                    StateList::addToList( nm1->content, l->current );
                      // remove State pointer from list to avoid deletion of
                      // moved State
                    l->current = NULL;
                      // go to next minimum merged with m2
                    l = l->next;
                }
            }
              // update tree structure
            Node* m1MergeParent = nm1;
            if (nm2->parent->content != NULL) {

                m1MergeParent 
                = getHighestParentLEQ( nm1 
                                , nm2->parent->content->current->getEnergy());
            } else {
                m1MergeParent = root;
            }
            
              // remove m2 from children list of m2->parent
            NodeList::remFromList( nm2->parent->children, nm2 );
            
              // consolidate changed tree
            Node* newRoot = mergeSubTrees( m1MergeParent, nm2->parent );
            
              // remove obsolete empty node nm2
            delete nm2;
              // remind the loss of a leaf
            numOfLeaves--;
            
              // check if new root has only one child and is therefore obsolete
            checkIfObsolete( newRoot );
            
              // barrier tree was changed due to minima (leave) merge
            return true;
        }

        
          // --> check if barrier is lower than already know one
        Node* knownBarrier = getLCA( nm1, nm2 );
        if (    knownBarrier->content == NULL  // in different trees of forest
                || State::less( s, knownBarrier->content->current ) 
            )
        {
            
              // check if known barrier has same energy as new state
              // --> node barrier state has to be updated
            if (knownBarrier->content != NULL
                && knownBarrier->content->current->getEnergy() == s->getEnergy())
            {
                  // because we are in this program part : s < knownBarrier
//                // remove old barrier state
//              delete knownBarrier->content->current;
                  // set new barrier state
                knownBarrier->content->current = s->clone(knownBarrier->content->current);
                  // barrier was changed
                return true;
            }

            // --> update barrier
              // get nodes to join below new barrier
            Node* t1 = getHighestParentLEQ( nm1, s->getEnergy());
            Node* t2 = getHighestParentLEQ( nm2, s->getEnergy());
              // get their parent nodes
            Node* pt1 = t1->parent;
            Node* pt2 = t2->parent;

              // pointer to newly generated node
            Node* newNode = NULL;
            
              // check if t1 represents already barrier height to add to tree
            if (t1->content->current->getEnergy() == s->getEnergy()) {
                  // no new node has to be added
                newNode = t1;
                  // exchange the barrier state if the new one is smaller
                  // according to the order
                if (State::less(s,newNode->content->current))
                {
//                    // remove old barrier state
//                  delete newNode->content->current;
                      // set new barrier state
                    newNode->content->current = s->clone(newNode->content->current);
                }
                
            } else {
                  // add new internal node containing s with children t1 and t2
                newNode = new Node( root
                                    , new NodeList(t1,NULL)
                                    , new StateList(s->clone(), NULL )
                                    );
            }

              // set new parent for t1 and t2
            NodeList::remFromList( pt1->children, t1 );
            NodeList::remFromList( pt2->children, t2 );
            t1->parent = newNode;
            t2->parent = newNode;
            
              // check second node if to join below barrier
            if (t2->content->current->getEnergy() 
                    < newNode->content->current->getEnergy()
                || t2->children == NULL  ) // is leaf 
            {
                  // we add the node to the children list
                NodeList::addToList( newNode->children, t2 );
            } else {
            
                 // t2 is internal node
                 // t2 is obsolete because it represents the same barrier height
                 // as newNode --> join both and remove t2
                
                  // exchange the barrier state if the new one is smaller
                  // according to the order
                if (State::less(    t2->content->current
                                    , newNode->content->current))
                {
                      // remove old barrier state
                    delete newNode->content->current;
                      // set new barrier state
                    newNode->content->current = t2->content->current;
                      // remove moved barrier state to avoid recursive deletion
                    t2->content->current = NULL;
                }
                    
                  // move all children hold by t2 into newNode
                NodeList* n = t2->children;
                while (n != NULL) {
                      // add all states that are merged with t2 to newNode
                    NodeList::addToList( newNode->children, n->current );
                      // update parent information
                    n->current->parent = newNode;
                      // remove State pointer from list to avoid deletion of moved
                      // State
                    n->current = NULL;
                      // go to next minimum merged with m2
                    n = n->next;
                }
                
                  // remove obsolete and empty t2
                delete t2;
            }
            
            
              // add new node as child to lower barrier from pt1/pt2
            if (pt1 == pt2) { 
                  // same sub tree root parent --> one operation sufficient
                NodeList::addToList( pt1->children, newNode );
                  // set parent node of newNode
                newNode->parent = pt1;
                  // check if pt1/pt2 is obsolete due to only a single child left
                pt1 = checkIfObsolete( pt1 );
            } else {
                if (    pt1 != root
                        && (pt2 == root 
                                || pt1->content->current->getEnergy() 
                                    < pt2->content->current->getEnergy() )
                    )
                {  // pt1 is the lower of both barriers
                      // add newNode to pt1 children
                    NodeList::addToList( pt1->children, newNode );
                      // set parent node of newNode
                    newNode->parent = pt1;
                      // merge the trees of pt1 and pt2

                    Node* newRoot = mergeSubTrees( pt1, pt2 );
                      // check if newRoot is obsolete due to only a single child left
                    newRoot = checkIfObsolete( newRoot );
                } else {
                    // pt2 is the lower of both barriers
                      // add newNode to pt2 children
                    NodeList::addToList( pt2->children, newNode );
                      // set parent node of newNode
                    newNode->parent = pt2;
                      // merge the trees of pt1 and pt2
                    Node* newRoot = mergeSubTrees( pt2, pt1 );
                      // check if newRoot is obsolete due to only a single child left
                    newRoot = checkIfObsolete( newRoot );
                }
            }

              // barrier was updated
            return true;
            
        }
        
          // known barrier is smaller than the new one to update 
        return false;
    }


    bool
    LT_BarrierTree::
	addSaddle(const State* const si, const State* const sj,const State* const s) 
    
    {
        
        bool hasChanged = false;
        
          // get index of the minimum i
        size_t idx_i = getMinIndex(si);
          // if not known add to tree
        if (idx_i == LandscapeTopology::INVALID_INDEX ) {
            idx_i = addMin(si);
            hasChanged = true;
        }
          // get index of the minimum j 
        size_t idx_j = getMinIndex(sj);
          // if not known add to tree
        if (idx_j == LandscapeTopology::INVALID_INDEX ) {
            idx_j = addMin(sj);
            hasChanged = true;
        }
        
          // add saddle and return if barrier tree was changed
        return addSaddle( idx_i, idx_j, s) || hasChanged;
    }

    const size_t 
    LT_BarrierTree::
	addMin(const State* const s) 
    {
        assertbiu(s != NULL, "given minimum State to add is NULL");
        assertbiu(getMinIndex(s) == LandscapeTopology::INVALID_INDEX
                , "given minimum is known and already part of the landscape");

          // make copy of the minimum to add
        State* newMin = s->clone();
        
          // add the minimum to the minima reference
        vMinima.push_back(newMin);
        
          // create new node to add to the tree
        Node* newNode = new Node( NULL, NULL, new StateList(newMin,NULL) );
          // update minima --> node mapping
        min2node[newMin] = newNode;
        
          // add node to the tree
        if (root == NULL) {
              // no tree available --> make new minimum to root
            root = newNode;
            newNode->parent = NULL;
        } else {
            if (root->content != NULL) {
                  // current root is a minimum or a barrier, but we have to
                  // create a forest with virtual root node
                Node* newRoot 
                        = new Node(NULL, new NodeList(root,NULL), NULL, true);
                root->parent = newRoot;
                root = newRoot;
            }
              // add new minimum as a new barrier tree to the forest
            NodeList::addToList(root->children, newNode);
            newNode->parent = root;
        }

         // update mfeIndex if necessary
        if (vMinima.size() == 1) {
            mfeIndex = 0;
        } else {
            if (State::less( s, vMinima[mfeIndex])) {
                  // new energy is smaller 
                  // OR new energy is equal but string representation is smaller
                mfeIndex = vMinima.size()-1;
                  // is not sorted anymore because smallest is at the end now
                sorted = false;
            } else if (sorted) {
                  // compare last with newly added if the new one is not smaller
                sorted = State::less( vMinima[vMinima.size()-2], s );
            }
        }
        
          // remind the newly added leaf
        numOfLeaves++;

        return vMinima.size()-1;
    }



    // calculates the distance between this and a given LT_BarrierTree
    double 
    LT_BarrierTree::
	getDistance(const LandscapeTopology* const pLT ) const 
    {
        assertbiu(pLT!=NULL, "given LandscapeTopology is NULL!");
        
        const LT_BarrierTree* const bt 
                    = dynamic_cast<const LT_BarrierTree* const>(pLT);

        if ( bt == NULL ) { 
            // no LT_SaddleNet --> use default function
            return LandscapeTopology::getDistance(pLT);
        }
            
        double squareDeviation = 0.0;
        size_t numOfDeviations = 0;
        double minimaDistance = 0.0;
        size_t numOfDifferentMinima = 0;
        
//      @TODO TESTING !!!
//for (std::map<const State*,Node*>::const_iterator it=min2node.begin();it!=min2node.end();it++) {
//  //std::cerr <<"## " <<getMinIndex(it->first) <<" : " <<it->first->toString() <<" --> " <<it->second <<" --> " <<it->second->parent <<" --> " <<it->second->parent->parent<<"\n";
//}
          // set up mapping of leaves (minima) between the two barrier trees
        Node2Node_MAP minMap1to2;
        Node2Node_MAP minMap2to1;
          // find and mark the intersection of both minima(leave) sets
        std::stack<const NodeList*> childStack;
          // find mapping for all leaves of this barrier tree
        if (this->root != NULL) {
            childStack.push(this->root->children);
        }
        const NodeList* nextChild = NULL;
        while( !childStack.empty() ) {
            nextChild = childStack.top();
            if (nextChild == NULL)  {
                  // last of the current children list was handled
                childStack.pop();
            } else {
                  // set stack top to the next in the current children list
                childStack.top() = nextChild->next;
                
                assertbiu( nextChild->current != NULL
                                , "child information has no node link");
                
                  // check children of current node
                if (nextChild->current->children != NULL) {
                      // internal node found
                    childStack.push( nextChild->current->children );
                } else {
                      // leaf found
                      // check if one of the minima of this leaf are present in
                      // the other barrier tree
                    StateList* s = nextChild->current->content;
                    while( s != NULL ) {
                          // check if the current minimum of the leaf is element
                          // of the other tree
                        size_t to = bt->getMinIndex(s->current);
                          // store leave/minima mapping
                        if (to != INVALID_INDEX) {
                              // get corresponding leaf
                            const Node* toNode = bt->min2node.find(
                                                    bt->vMinima[to])->second;
                              // update mapping
                            minMap1to2[nextChild->current] = toNode;
                            minMap2to1[toNode] = nextChild->current;
                              // found mapping --> leave minima loop
                            break;
                        }
                          // go to next minimum of this leaf
                        s = s->next;
                    }
                    if (minMap1to2.find(nextChild->current) != minMap1to2.end())
                    { // still not present in m2 --> handle missing minimum 
                        numOfDifferentMinima++;
                        minimaDistance += MISSING_MINIMUM_PENALTY 
                                * exp( -(nextChild->current->content
                                            ->current->getEnergy())
                                        / BOLTZMANN_KT );
                    }
                }
            }
        }
          // find mapping for all leaves of the other barrier tree
        if (bt->root != NULL) {
            childStack.push(bt->root->children);
        }
        while( !childStack.empty() ) {
            nextChild = childStack.top();
            if (nextChild == NULL)  {
                  // last of the current children list was handled
                childStack.pop();
            } else {
                  // set stack top to the next in the current children list
                childStack.top() = nextChild->next;
                
                  // check children of current node
                if (nextChild->current->children != NULL) {
                      // internal node found
                    childStack.push( nextChild->current->children );
                } else {
                      // leaf found --> check if already mapped
                    if (minMap2to1.find(nextChild->current) != minMap2to1.end()) 
                    { 
                          // check if one of the minima of this leaf are present
                          // in the other barrier tree
                        StateList* s = nextChild->current->content;
                        while( s != NULL ) 
                        {
                              // check if the current minimum of the leaf is element
                              // of the other tree
                            size_t to = this->getMinIndex(s->current);
                              // store leave/minima mapping
                            if (to != INVALID_INDEX) {
                                  // get corresponding leaf
                                const Node* toNode = this->min2node.find(
                                                    this->vMinima[to])->second;
                                  // update mapping
                                minMap2to1[nextChild->current] = toNode;
                                  // found mapping --> leave minima loop
                                break;
                            }
                              // go to next minimum of this leaf
                            s = s->next;
                        }
                    }
                    if (minMap2to1.find(nextChild->current) != minMap2to1.end()) 
                    { // still not present in m1 --> handle missing leaf 
                        numOfDifferentMinima++;
                        minimaDistance += MISSING_MINIMUM_PENALTY 
                                * exp( -(nextChild->current->content
                                            ->current->getEnergy())
                                        / BOLTZMANN_KT );
                    }
                }
            }
        }
        
          // normalize minima distance
        minimaDistance /= (double)numOfDifferentMinima;

        // calculate barrier square deviation
        
        // --> get maximal barrier height maxDelta of this and the other tree
        // this value will be used if a barrier in one of the trees is missing
        double maxBarrierHeight = 0.0;
        {
            // NOTE: only ensuringing highest barrier if this is no forest ...
            
              // get max barrier height of this tree
            Node* mfeRoot = this->min2node.find(this->getMFEState())->second;
            while( mfeRoot->parent != NULL && mfeRoot->parent->content != NULL ) {
                mfeRoot = mfeRoot->parent;
            }
            if (mfeRoot->content != NULL) { // mfe is no single tree in forest
                maxBarrierHeight = mfeRoot->content->current->getEnergy()
                                    - this->getMFEState()->getEnergy();
            }
              // get max barrier height of the other tree
            mfeRoot = bt->min2node.find(bt->getMFEState())->second;
            while( mfeRoot->parent != NULL && mfeRoot->parent->content != NULL ) {
                mfeRoot = mfeRoot->parent;
            }
            if (mfeRoot->content != NULL) { // mfe is no single tree in forest
                  // check if this barrier height is higher than of this tree
                maxBarrierHeight = std::max( maxBarrierHeight 
                                        , mfeRoot->content->current->getEnergy()
                                            - bt->getMFEState()->getEnergy());
            }
        }
          // calculate square value
        maxBarrierHeight *= maxBarrierHeight;
        
        // --> calculate square deviations of all barriers between minima leaves
        //     that are present in both trees
          // needed to know if a minimum leaf is unique in the second tree
        std::set<const Node*> mappedFrom1;
          // check all barriers present in this tree
        typedef Node2Node_MAP::const_iterator MapIter;
        for (MapIter from=minMap1to2.begin(); from != minMap1to2.end(); from++) 
        {
              // store that this node was reached from this tree
            mappedFrom1.insert(from->second);
              // calc barriers to all other leaves
            for (MapIter to = from; ++to != minMap1to2.end(); ) 
            {
                  // calculate barrier
                const Node* lca1 = this->getLCA(from->first, to->first);
                const Node* lca2 = bt->getLCA(from->second, to->second);
                if (lca1->content != NULL && lca2->content != NULL) {
                    // there exists a barrier in both trees 
                    // --> calculate square deviation
                      // get barrier height of 'from' node in this tree 
                    double bDiff = from->first->content->current->getEnergy()
                                    - lca1->content->current->getEnergy();
                      // substract barrier height of 'from' node in the other tree 
                    bDiff -= from->second->content->current->getEnergy()
                                - lca2->content->current->getEnergy();
                      // calculate square deviation and store
                    squareDeviation += (bDiff * bDiff);
                } else
                if (    (lca1->content != NULL && lca2->content == NULL)
                        || (lca1->content == NULL && lca2->content != NULL) ) 
                {
                    // a MISSING BARRIER found --> increase deviation
                    squareDeviation += maxBarrierHeight;
//              } else {
                  // both are in different trees in the forest 
                  // --> therefore no barrier difference and no squareDeviation
                }
                  // count that we have handled a barrier
                numOfDeviations++;
            }
        }
          // check all remaining barriers present in the other tree ONLY
        for (MapIter from=minMap2to1.begin(); from != minMap2to1.end(); from++) 
        {
              // check if this leaf was already mapped
            if (mappedFrom1.find(from->first) != mappedFrom1.end()) {
                continue;
            }
              // calc barriers to all other leaves
            for (MapIter to = minMap2to1.begin(); to != minMap2to1.end(); to++) 
            {
                  // check if this is the source and thus to skip
                if (to->first == from->first) {
                    continue;
                }
                
                  // calculate barrier
                const Node* lca2 = bt->getLCA(from->first, to->first);
                const Node* lca1 = this->getLCA(from->second, to->second);
                if (lca1->content != NULL && lca2->content != NULL) {
                    // there exists a barrier in both trees 
                    // --> calculate square deviation
                      // get barrier height of 'from' node in this tree 
                    double bDiff = from->first->content->current->getEnergy()
                                    - lca1->content->current->getEnergy();
                      // substract barrier height of 'from' node in the other tree 
                    bDiff -= from->second->content->current->getEnergy()
                                - lca2->content->current->getEnergy();
                      // calculate square deviation and store
                    squareDeviation += (bDiff * bDiff);
                } else
                if (    (lca1->content != NULL && lca2->content == NULL)
                        || (lca1->content == NULL && lca2->content != NULL)) 
                {
                    // a MISSING BARRIER found --> increase deviation
                    squareDeviation += maxBarrierHeight;
//              } else {
                  // both are in different trees in the forest
                  // --> therefore no barrier difference and no squareDeviation
                }
                  // count that we have added a barrier
                numOfDeviations++;
            }
        }
        
        
          // nothing was matched etc --> use minima distance only
        if (numOfDeviations == 0) {
            return minimaDistance;
        }
          // calculate distance via barrier RMSD and minima distance
        return minimaDistance
                + sqrt( squareDeviation / (double)numOfDeviations );
    }


    std::pair<int,std::string>
    LT_BarrierTree::read(   std::istream & in, 
                        const State * const templateState )
    {

        typedef std::pair<int,std::string> RETURNTYPE;
        State* s;
        
          // stores the mapping of the minima indices in the file to the indices
          // used in the barrier tree to add the barriers
        Int2SizeT_MAP idx2min;
        
        std::string line;
        std::stringstream sstr;
        
        bool headerParsed = false;
        size_t dim = 0;
        size_t lineNumber = 0;
        
        std::string::size_type i;
        

        while( std::getline(in, line) ) {
            ++lineNumber;
            s = NULL;
            
            if(line.find("GRAPHVIZ") != std::string::npos) {
                break;
            }
            
              // check if next line describes the header information
            i = line.find(OUTPUT_KEY_ELL_HEAD);
            if(i != std::string::npos) {
                
                std::string lt_id("");
                dim = 0;
                std::string stateDescr("");
                
                  // parse the minimum
                RETURNTYPE ret 
                    = LandscapeTopology::parseHeader(   line.substr(i), 
                                                        lt_id,
                                                        dim,
                                                        stateDescr );
                  // check for parsing errors
                if (ret.first != 0) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : "<<ret.second;
                    return RETURNTYPE(ret.first,sstr.str());
                }
                  // check for semantical errors
                if ( lt_id.compare(LT_MinimaSet::LT_ID) != 0 
                     && lt_id.compare(LT_BarrierTree::LT_ID) != 0
                     && lt_id.compare(LT_SaddleNet::LT_ID) != 0) 
                {
                    return RETURNTYPE(-1,"ELL-LT-TYPE not supported");
                }
                if (stateDescr.compare(templateState->getID()) != 0 ) {
                    return RETURNTYPE(-2,"STATETYPE differs from given type");
                } 
                headerParsed = true;
            }
            
              // check if next line describes a minimum
            i = line.find(OUTPUT_KEY_MINIMUM);
            if(i != std::string::npos) {
                  // check if header already parsed
                if (!headerParsed) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : minimum present but no header found so far";
                    return RETURNTYPE(-2,sstr.str());
                }
                
                  // count the new found minimum
                if (dim == 0) {
                    return RETURNTYPE(-3,"MORE minima present than given via DIMENSION");
                }
                dim--;
                
                  // parse the minimum
                RETURNTYPE ret 
                    = LandscapeTopology::parseMinimum(  *this, 
                                                        line.substr(i), 
                                                        templateState, 
                                                        idx2min );
                
                  // check for parsing errors
                if (ret.first != 0) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : "<<ret.second;
                    return RETURNTYPE(ret.first,sstr.str());
                }
                  // check for minima merged with the added one
                if (ret.second.size() > 0) {
                      // get leaf node of last minimum added
                    Node* lastAddedNode = min2node[vMinima[vMinima.size()-1]];
                      // parse remaining minima
                    size_t start = ret.second.find_first_not_of(' ',0);
                    while (start < ret.second.size()) {
                          // get end of encoding
                        size_t end = ret.second.find_first_of(' ',start);
                        State* subMin = templateState->fromString(ret.second.substr(start,end-start));
                        assertbiu(subMin != NULL, "minimum parsing failed");
                          // add min to minima list of last added node
                        StateList::addToList( lastAddedNode->content, subMin);
                          // add min to global minima list
                        vMinima.push_back(subMin);
                        min2node[subMin] = lastAddedNode;
                          // go to next minimum encoding
                        start = ret.second.find_first_not_of(' ',end);
                    }
                }
            } 
            
              // check if next line describes a barrier
            i = line.find(OUTPUT_KEY_BARRIER);
            if(i != std::string::npos) {
                  // check if header already parsed
                if (!headerParsed) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : barrier present but no header found so far";
                    return RETURNTYPE(-2,sstr.str());
                }
                
                  // parse the barrier
                RETURNTYPE ret 
                    = LandscapeTopology::parseBarrier(  *this, 
                                                        line.substr(i), 
                                                        templateState, 
                                                        idx2min );
                
                  // check for parsing errors
                if (ret.first != 0) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : "<<ret.second;
                    return RETURNTYPE(ret.first,sstr.str());
                }
            } 
            
              // check if next line describes a saddle
            i = line.find(OUTPUT_KEY_SADDLE);
            if(i != std::string::npos) {
                  // check if header already parsed
                if (!headerParsed) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : saddle present but no header found so far";
                    return RETURNTYPE(-2,sstr.str());
                }
                
                  // parse the saddle
                RETURNTYPE ret 
                    = LandscapeTopology::parseSaddle(   *this, 
                                                        line.substr(i), 
                                                        templateState, 
                                                        idx2min );
                
                  // check for parsing errors
                if (ret.first != 0) {
                    sstr.clear();
                    sstr <<"line " <<lineNumber
                        <<" : "<<ret.second;
                    return RETURNTYPE(ret.first,sstr.str());
                }
            } 

              // garbage collection
            if (s != NULL) {
                delete s;
                s = NULL;
            }
            
        }

        if (dim > 0) {
            return RETURNTYPE(-3,"LESS minima present than given via DIMENSION");
        }
        
        return RETURNTYPE(0,"");
    }
    
    void
    LT_BarrierTree::
	write( std::ostream& out, 
            const bool writeGraph ) const
    {
        const size_t dim = numOfLeaves;


        const std::string dummyID = "ell::State";
        LandscapeTopology::writeHeader(out, LT_ID, dim, (vMinima.size()==0?dummyID:vMinima[0]->getID()));


          // print barriers
        Node2Node_MAP node2firstLeaf;
        Node2SizeT_MAP node2idx;
        SizeT2Node_MAP idx2node;
        size_t internalNodeCount = dim-1;
        size_t leafNodeCount = 0;
        for (size_t i=0; i<vMinima.size(); i++) {
            Node* cur = min2node.find(vMinima[i])->second;
              // check if leaf was already handled
            if (node2idx.find(cur) != node2idx.end()) {
                continue;
            }
              // label node
            node2idx[cur] = leafNodeCount;
            idx2node[leafNodeCount] = cur;
              // print smallest minimum stored in this leaf
            if (cur->content->next == NULL) {
                LandscapeTopology::writeMinimum( out, cur->content->current, leafNodeCount );
            } else {
                std::vector<const State*> minima;
                const StateList* next = cur->content;
                while (next != NULL) {
                    minima.push_back(next->current);
                    next = next->next;
                }
                LandscapeTopology::writeMinimum( out, minima, leafNodeCount );
            }
              // increase node counter
            leafNodeCount++;
              // prepare parent recursion
            node2firstLeaf[cur] = cur;
            bool nextRecursion = true;
            while( nextRecursion && cur->parent != NULL ) {
                if (node2idx.find(cur->parent)==node2idx.end()) {
                    internalNodeCount++;
                    node2idx[cur->parent] = internalNodeCount;
                    idx2node[internalNodeCount] = cur->parent;
                    node2firstLeaf[cur->parent] = node2firstLeaf[cur];
                } else {
                      // have found a known internal node 
                      // --> no recursion needed
                    nextRecursion = false;
                }
                  // go recursive a level up in the tree
                cur = cur->parent;
            }
        }
        assertbiu(numOfLeaves==leafNodeCount, "not all leaves printed");
          // vector that holds the minima connected by a barrier
        std::vector<size_t> minConnected;
          // write barriers to stream
        for (SizeT2Node_MAP::iterator it = idx2node.begin();
            it != idx2node.end(); it++)
        {

              // check if internal node and no virtual root of a forest
            if (it->second->children != NULL && it->second->content != NULL) {
                  // reset minima collector
                minConnected.clear();
                  // collect connected minima
                NodeList* children = it->second->children;
                while (children != NULL) {
                    minConnected.push_back(
                            node2idx[ node2firstLeaf[ children->current ] ] );
                    children = children->next;
                }
                
                  // write barrier to stream
                LandscapeTopology::
					writeBarrier(  out,
                                    it->second->content->current,
                                    it->first,
                                    minConnected );
            }
        }
    
        if (writeGraph) {
        
            out << "\n" << "//GRAPHVIZ"<<"\n";
            
            out << "digraph BARRIERTREE {"<<"\n";
            for (size_t i = 0; i < dim; ++i) {
                out << "\t\tM"<<i<<";\n";
            }
            
              // print internal barrier nodes and edges
            for (SizeT2Node_MAP::iterator it = idx2node.begin();
                it != idx2node.end(); it++)
            {
                if (it->second == root) {
                    out << "\t\tB"<<node2idx[root]<<"; /* root */\n";
                }
                  // check if internal node and no virtual root of a forest
                if (it->second->children != NULL && it->second->content != NULL) {
                      // print barrier index
                    // print connected minima
                    NodeList* children = it->second->children;
                    while (children != NULL) {
                        size_t curIdx = node2idx[children->current];
                        out << "\t\tB"<<it->first<<" -> "
                            <<(curIdx<dim?"M":"B")<<curIdx<<";\n";
                        children = children->next;
                    }
                }
            }
            
            out << "label = \"\\n\\nLT_BarrierTree Diagram\\nwith "<<dim
                <<" minima \";\n"
                << "fontsize=20;\n"
                << "}\n";
        }
        
          // flush the output stream
        out.flush();
        
    }
    
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	getHighestParentLEQ( Node* child, const double maxEnergy )
    {
        assertbiu( child != NULL, "starting node is NULL" );
        
        Node* ret = child;
        
          // find (grand)parent node of child that is either root or internal 
          // node with energy below or equal to threshold 
        while(  ret != root
                && ret->parent->content != NULL
                && ret->parent->content->current->getEnergy() <= maxEnergy ) 
        {
            ret = ret->parent;
        }
        
        return ret;
    }
    
    const 
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	getHighestParentLEQ( const Node* child, const double maxEnergy ) const
    {
        assertbiu( child != NULL, "starting node is NULL" );
        
        const Node* ret = child;
        
          // find (grand)parent node of child that is either root or internal 
          // node with energy below or equal to threshold 
        while(  ret != root
                && ret->parent->content != NULL
                && ret->parent->content->current->getEnergy() <= maxEnergy ) 
        {
            ret = ret->parent;
        }
        
        return ret;
    }
    
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	mergeSubTrees( Node* ta, Node* tb )
    {
        if (ta == tb) {
            return ta;
        }
          // check if merge is needed
        if (    ta->content != NULL
                && tb->content != NULL
                && ta->content->current->getEnergy() 
                    == tb->content->current->getEnergy() )
        {
              // t1 has lower energy than t2
            Node* t1 = ta;
            Node* t2 = tb;
            // same barrier height --> merge nodes
              // --> ensure that (t1->current < t2->current)
            if (State::less(t2->content->current, t1->content->current)) {
                t1 = tb;
                t2 = ta;
            }
              // --> merge t2 into t1
            assertbiu(t2->children != NULL, "t2 is no internal node");
              // list of nodes that are children of t2 --> move to t1
            {
                NodeList*& n = t2->children;
                while( n != NULL ) {
                      // add all minima that are merged with t2 to t1
                    NodeList::addToList( t1->children, n->current );
                      // update parent information of new children
                    n->current->parent = t1;
                      // remove node pointer from list to avoid deletion of
                      // moved nodes
                    n->current = NULL;
                      // go to next minimum merged with t2
                    n = n->next;
                }
            }
            
              // --> remove t2 from t2->parent->children
            NodeList::remFromList( t2->parent->children, t2 );
              // store t2 parent information
            Node* pt2 = t2->parent;
              // --> delete t2
            delete t2;
              // --> return merge of t1->parent an t2->parent
            return mergeSubTrees(t1->parent,pt2);
        }

          // check if nodes are not obsolete
        if (ta->children != NULL) {
            ta = checkIfObsolete(ta);
        }
        if (tb->children != NULL) {
            tb = checkIfObsolete(tb);
        }
        
          // check if no more merge neccessary
        if (ta == tb || ta == tb->parent || ta == root) {
            return ta;
        } 
        if (tb == ta->parent || tb == root) {
            return tb;
        }
        
          // t1 has lower energy than t2
        Node* t1 = ta;
        Node* t2 = tb;
        if (    t1->content->current->getEnergy() 
                > t2->content->current->getEnergy() ) 
        {
            t1 = tb;
            t2 = ta;
        }
        
          // store old parents
        Node* pt1 = t1->parent;
        Node* pt2 = t2->parent;
        
          // check if t1 has a parent that is still below t2
        if (    pt1->content != NULL 
                && State::less(pt1->content->current, t2->content->current)) 
        {
              // merge parent of t1 with t2
            return mergeSubTrees( pt1, t2);
        }
        
          // --> insert t2 between t1 and t2->parent
        NodeList::remFromList( pt1->children, t1 );
        NodeList::addToList( t2->children, t1 );
        t1->parent = t2;
        NodeList::addToList( pt1->children, t2 );
        t2->parent = pt1;
          // --> remove t2 from t2->parent->children
        NodeList::remFromList( pt2->children, t2 );
        
          // recursive call with new joint subtrees
        return mergeSubTrees( pt2, pt1 );
    }
    

    const 
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	getLCA( const Node* n1, const Node* n2 ) const 
    {
          // simple checks for LCA
        if (n1 == n2 )
            return n1;
        if (n1->parent == n2->parent)
            return n1->parent;
        
        const Node* lca = root;
        
          // path from n1/n2 to root node
        std::list<const Node*> p1, p2;
        
          // find path to root from n1
        p1.push_back(n1);
        while( (*(p1.rbegin())) != root ) {
            p1.push_back((*(p1.rbegin()))->parent);
        }
          // find path to root from n2
        p2.push_back(n2);
        while( (*(p2.rbegin())) != root ) {
            p2.push_back((*(p2.rbegin()))->parent);
        }
        
        std::list<const Node*>::const_reverse_iterator i1 = p1.rbegin();
        std::list<const Node*>::const_reverse_iterator i2 = p2.rbegin();
        
        while(*i1 == *i2) {
            lca = *i1;
            i1++;
            i2++;
        }
        
        return lca;
    }
    
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	getLCA( Node* n1, Node* n2 ) 
    
    {
          // simple checks for LCA
        if (n1 == n2 )
            return n1;
        if (n1->parent == n2->parent)
            return n1->parent;
        Node* lca = root;
        
          // path from n1/n2 to root node
        std::list<Node*> p1, p2;
        
          // find path to root from n1
        p1.push_back(n1);
        while( (*(p1.rbegin())) != root ) {
            p1.push_back((*(p1.rbegin()))->parent);
        }
          // find path to root from n2
        p2.push_back(n2);
        while( (*(p2.rbegin())) != root ) {
            p2.push_back((*(p2.rbegin()))->parent);
        }
        
        std::list<Node*>::const_reverse_iterator i1 = p1.rbegin();
        std::list<Node*>::const_reverse_iterator i2 = p2.rbegin();
        
        while(*i1 == *i2) {
            lca = *i1;
            i1++;
            i2++;
        }
        
        return lca;
    }
    
    
    const State*
    LT_BarrierTree::
	getBarrier( const size_t m1, const size_t m2 ) const 
    {
        assertbiu(m1 < vMinima.size(), "m1 index out of range (vMinima)");
        assertbiu(m2 < vMinima.size(), "m2 index out of range (vMinima)");
        
          // get LCA node
        const Node* lca = getLCA(   min2node.find(vMinima[m1])->second
                                    , min2node.find(vMinima[m2])->second );
        
          // get barrier state or NULL if in non-connected trees of a forest
        return lca->content->current;
    }
    
    
    const State*
    LT_BarrierTree::
	getBarrier( const State* const m1, const State* const m2 ) const 
    {
        const size_t im1 = getMinIndex(m1);
        const size_t im2 = getMinIndex(m2);
        
        assertbiu(im1 != INVALID_INDEX, "m1 is not part of the topology");
        assertbiu(im2 != INVALID_INDEX, "m2 is not part of the topology");

          // calculate barrier
        return getBarrier( im1, im2 );
    }
    
    LT_BarrierTree::
    Node*
    LT_BarrierTree::
	checkIfObsolete( Node* n )
    {
//std::cerr <<"###### CHECK IF OBSOLETE ( " <<n <<" )";
        assertbiu( n->children != NULL, "node n has no children (NULL)");
        assertbiu( n->children->current != NULL, "node n has no valid children (current==NULL)");

          // return value
        Node* finalRoot = n;
        
          // check if only one child present
        if (n->children->next == NULL) {
            if (n == root) {
                  // set root a level lower
                root = n->children->current;
                  // update parent information of the new root node
                root->parent = NULL;
                  // set final root
                finalRoot =  root;
            } else {
                  // remove n from parents children list
                NodeList::remFromList(n->parent->children, n);
                  // add the only child of n to parents children list
                NodeList::addToList(    n->parent->children
                                        , n->children->current);
                  // update parent of child of n
                n->children->current->parent = n->parent;
                  // set final root
                finalRoot = n->parent;
            }
              // erase connection of n to tree
            n->children->current = NULL;
              // delete newRoot
            delete n;
//std::cerr <<" = TRUE\n";
        } else {
//std::cerr <<" = false\n";
        }
        
        return finalRoot;
    }
    

    bool
    LT_BarrierTree::
	isForest( void ) const
    {
        return (root!=NULL) && (root->content == NULL);
    }

    

    size_t
    LT_BarrierTree::
	getNumOfTrees( void ) const
    {
        size_t trees = 0;
        
          // check if forest
        if (root->content == NULL) {
            NodeList* child = root->children;
            while (child != NULL) {
                trees++;
                child = child->next;
            }
        } else {
            trees = 1;
        }
        
        return trees;
        return (root!=NULL) && (root->content == NULL);
    }

    
    size_t
    LT_BarrierTree::
	getNumOfLeaves( void ) const
    {
        assertbiu(numOfLeaves <= vMinima.size(), "there is an inconsistency : more leaves than minima present");
        
        return numOfLeaves;
    }

    
    void 
    LT_BarrierTree::
	write_PStree( std::ostream& out, const size_t maxLeaves2print) 
    
    {
        if(this->isForest()) {
#ifdef DEBUG
            std::cerr << "SORRY : this is a forest which is currently not supported for PStree .. ";
#endif
            return;
        }
        
          // sort minima
        this->sort();
          // get maximal number of leaves to print
        const size_t max_print = std::min( maxLeaves2print, numOfLeaves);
        
        std::vector<LT_BarrierTree::nodeT> nodes(max_print);

          // all leaves are attracted by the mfe of the tree
        
        Node2Int_MAP node2idx;

        int curMin = 0;
        for( int addedMin=0; addedMin < (int)max_print; addedMin++ ) {
            
              // go to next minimum that represents a leaf
            while (min2node[vMinima[curMin]]->content->current != vMinima[curMin]) {
                curMin++;
            }

            {
                Node* node = min2node[vMinima[curMin]];
                while (node != NULL && node2idx.find(node)==node2idx.end()) {
                    node2idx[node] = addedMin;
                    node = node->parent;
                }
            }
            
            register int f;
            
            Node* parent = min2node[vMinima[curMin]]->parent;
            assertbiu(parent->content != NULL && parent->content->current != NULL, "is a forest");
            double E_saddle = parent->content->current->getEnergy();
            
            assertbiu(node2idx.find(parent) != node2idx.end(), "father not found");
            f = node2idx[parent];

            nodes[addedMin].father = (f==addedMin)?-1:f;
            
            nodes[addedMin].height = vMinima[curMin]->getEnergy();
            nodes[addedMin].saddle_height = E_saddle;
            
//          if (print_labels) {
//              char *L;
//              char *s;
//              s = unpack_my_structure(Lmin[ii].structure);
//              if ((POV_size)&&(Lmin[ii].global)) {
//                  L = (char *) space(sizeof(char)*(3+strlen(s)));
//                  strcat(L,s); strcat(L," *");
//                  nodes[s1-1].label = L;
//              }
//              else
//                  nodes[s1-1].label = strdup(s);
//              free(s);
//          }
            nodes[addedMin].label = NULL;
            curMin++;
        }
        
//std::cerr<<"\n max_print = "<<max_print <<"\n";       
//for (size_t x=0; x<max_print; x++) {
//  std::cerr <<" node " <<x <<" ( h=" ;
//  std::cerr <<nodes[x].height; 
//  std::cerr <<" s=" ;
//  std::cerr <<nodes[x].saddle_height; 
//  std::cerr <<" f="; 
//  std::cerr <<nodes[x].father <<" )\n";
//}
//std::cerr<<"\n";      
        
        PS_tree_plot( out, nodes);
    }



    void
    LT_BarrierTree::
    PS_tree_plot(   std::ostream& out
                    , const std::vector<nodeT> &nodes
                    ) 
    {

        /* plot tree with leaf-nodes nodes, to file filename (stdout if NULL) */
        int i, k, ll, f=0;
        linkT *chain, *l;
        int bbox[4] = {72, 144, 522, 700};
        const int n = (int)nodes.size();

        
        /* make index list, sorted by saddle height */
        std::vector<int> sindex(n,0);
        for (i=0; i<n; i++) sindex[i]=i;
        cmp_saddle compare(&nodes);
        std::sort(sindex.begin(), sindex.end(), compare);
        
        /* make order (x-coordinates) of leaves, we use a linked list for this */
        chain = new linkT[n];
        /* start connecting the links of the chain */
        for (i=0; i<n; i++) { /* n-1 merges */
            k = sindex[i]; 
            f = nodes.at(k).father; /* ith saddle merges k and f */
            if (f==-1) f=0;
            if (k==f) continue;  /* lowest node doesn't merge */
            for (l= &chain[f]; l->next!=NULL; l = l->next); 
            l->next=&chain[k]; /* attach child to chain of father */
        }
        /* chain[f] now starts the ordered chain, next fill in the num field */
        for (i=0, l= &chain[f]; l!=NULL; l = l->next) 
            l->x = i++;
        
          // get current time
        time_t rawtime;
        struct tm * timeinfo;
        time ( &rawtime );
        timeinfo = localtime ( &rawtime );
        
        out <<
        "%!PS-Adobe-2.0 EPSF-1.2\n"
        "%Title: TreePlot\n"
        "%Creator: ell::LT_BarrierTree.cc\n"
        "%CreationDate: "
        <<asctime(timeinfo);
        out<<"%%BoundingBox: "
        <<bbox[0]<<" " << bbox[1] <<" " << bbox[2] <<" " << bbox[3] <<"\n";
        out <<
        "%%EndComments\n"
        "%%BeginProlog\n"
        "/treedict 100 dict def\n"
        "treedict begin\n"
        "% x y  => min(x,y)\n"
        "  /min { 2 copy gt { exch } if pop } bind def\n"
        "  /max { 2 copy lt { exch } if pop } bind def\n"
        "  /cmtx matrix currentmatrix def\n"
        "  /STR 128 string def\n"
        "  /NumH 1 def\n"
        "% - => -\n"
        "  /Init {\n"
        "    /LX [\n"
        "      LEAF {0 get} forall\n"
        "    ] def\n\n"
        "    /Helvetica findfont fsize scalefont setfont\n"
        "    /Lo [\n"
        "      (X) stringwidth pop % width\n"
        "      newpath 0 0 moveto\n"
        "      (X) true charpath\n"
        "      flattenpath pathbbox\n"
        "      pop exch pop exch sub neg 2 div % height\n"
        "     ] def\n"
        "  } def\n"
        "% - => -\n"
        "  /DrawScale {\n"
        "  gsave \n"
        "    maxy miny sub 30 div dup maxy add /maxy exch def miny sub /miny def\n"
        "    maxy miny sub log 0.9 sub floor 10 exch exp /tick exch def\n"
        "    newpath\n"
        "    LEAF length 0.5 sub 0 translate 0 miny moveto 0 maxy miny sub rlineto\n"
        "    miny tick div ceiling tick mul dup 0 exch moveto \n"
        "    maxy exch sub tick div cvi 1 add dup { % draw minor ticks\n"
        "      0.15 0 rlineto\n"
        "      -0.15 tick rmoveto\n"
        "    } repeat\n"
        "    % calculate major tick spacing (10, 5, or 2 minor ticks)\n"
        "    dup 69 gt { pop 10\n"
        "    } {\n"
        "      32 gt { 5 }\n"
        "      {2} ifelse\n"
        "    } ifelse\n"
        "    tick mul /mtick exch def\n"
        "    miny mtick div ceiling mtick mul dup 0 exch moveto\n"
        "    maxy exch sub mtick div cvi 1 add {\n"
        "      0.3 0 rlineto \n"
        "      gsave currentpoint 10 mul round 10 div cmtx setmatrix\n"
        "      STR cvs dup stringwidth pop\n"
        "      Lo aload pop 3 1 roll add neg exch rmoveto show pop\n"
        "      grestore\n"
        "      -0.3 mtick rmoveto\n"
        "    } repeat\n"
        "    cmtx setmatrix stroke    \n"
        "  grestore\n"
        "  } def\n"
        "% - => -\n"
        "  /SetBarFont {\n"
        "    matrix currentmatrix cmtx setmatrix\n"
        "    /Helvetica findfont fbsize scalefont setfont\n"
        "    setmatrix\n"
        "  } bind def\n"
        "% - => -\n"
        "  /SetLabelFont {\n"
        "    matrix currentmatrix cmtx setmatrix\n"
        "    /Courier findfont fsize scalefont setfont\n"
        "    setmatrix\n"
        "  } bind def\n"
        "% str => -\n"
        "  /Rotshow {\n"
        "    gsave\n"
        "      cmtx setmatrix -90 rotate\n"
        "      Lo aload pop\n"
        "      rmoveto show\n"
        "    grestore\n"
        "  } def\n"
        "% dy => - \n"
        "  /Rlineto {\n"
        "    dup abs MinHeight ge { % draw height at middle of line\n"
        "      dup gsave\n"
        "   dup 2 div 0 exch rmoveto\n"
        "   cmtx setmatrix -90 rotate\n"
        "   abs STR cvs dup stringwidth pop 2 div neg\n"
        "   //NumH rmoveto\n"
        "   show\n"
        "      grestore\n"
        "    } if\n"
        "    0 exch rlineto\n"
        "  } def\n"
        "% - => -\n"
        "  /Drawlabels {\n"
        "   0 LEAF {\n"
        "      aload pop moveto\n"
        "      dup LABEL exch get STR cvs Rotshow\n"
        "      1 add\n"
        "    } forall pop\n"
        "  } def\n"
        "% n => n'    Detect whether a minimum is connected\n"
        "  /MRX {\n"
        "     /murxi { true } def\n"
        "     dup 0 lt { pop 0 /murxi { false } def } if\n"
        "  } def\n"    
        "% - => -\n"
        "  /Connectlmins {\n"
        "    newpath\n"
        "    SADDEL {\n"
        "      /forest {false} def  %  draw as tree or forest node\n"
        "      aload pop exch dup 0 lt { pop 0 /forest {true} def} if"
        "   % => c h f\n"
        "      dup LX exch get [ exch LX 5 index get add 2 div "
        "% => c h f [ nx\n"
        "      3 index ]\t\t\t\t         % => c h f [ nx h ]\n"
        "      3 -1 roll dup LEAF 6 -1 roll get aload pop "
        "% => f [nx h] h h cx cy\n"
        "      dup 3 1 roll moveto\t\t         % => f [] h h cy\n"
        "      sub Rlineto                                % => f [] h\n"
        "      LEAF 3 index get aload pop exch\t\t % => f [] h fy fx\n"
        "      2 index forest {moveto} {lineto} ifelse \n"
        "      sub neg Rlineto\t\t\t         % => f [] h fy\n"
        "      LEAF 3 1 roll put\n"
        "    } forall\n"
        "    gsave\n"
        "      cmtx setmatrix stroke\n"
        "    grestore\n"
        "  } def\n"
        "% data starts here!!!\n"
        "  /LABEL [";

        /* print label array */
        if(nodes.at(0).label==NULL) 
            ll = 8;
        else 
            ll = 3+strlen(nodes.at(0).label);
        if(ll<8) 
            ll = 8;
        ll = (int) (80/ll);
        for (i=0; i<n; i++) {
            if (i%ll == 0)  
                out << "\n   ";
            if (nodes.at(i).label) 
                out << "(" << nodes.at(i).label <<") ";
            else
                out <<std::setw(3) <<(i+1) <<" ";
        }
        out <<"\n  ] def\n";

        /* print leaf node coordinates */
        out <<"% leaf node coordinates\n"
            <<"  /LEAF [";
        for (i=0; i<n; i++) {
            if (i%5 == 0)  out<< "\n   ";
            out<< "["<<std::setw(3)<<(chain[i].x)<<" "<<std::setw(10)<<std::fixed<<std::setprecision(3)<<nodes.at(i).height<<"] ";
        }
        out<< "  \n] def\n";

        /* print internal node coordinates */
        out<< "% internal nodes (saddle) coordinates, sorted by height\n"
        "  /SADDEL [";
        for (i=0; i<n; i++) {
            if (i%4 == 0)  out<< "\n   ";
            k=sindex[i]; 
            if (k==nodes.at(k).father) 
                continue;
            int fath = std::max(nodes.at(k).father,0);
            out<< "["<<std::setw(3)<<k <<" " <<std::setw(3)<<fath <<" "<<std::setw(10)<<std::fixed<<std::setprecision(3)<<nodes.at(k).saddle_height <<"] ";
        }
        delete chain;
        out <<
        "  \n] def\n"
        "end\n";
        out <<
        "%%EndProlog\n"
        "treedict begin\n"
        "  /fsize 10 def\n"
        "  /fbsize 7 def\n"
        "  Init\n"
        "  "<<(bbox[2]-1)<<" "<<bbox[1]<<" fsize 1.5 mul add translate\n";
        out<< "  "<<bbox[0]<<" "<<(bbox[2]-1)<<" sub LEAF length div % x-scale\n";
        out<< "  "<<(bbox[3]-1)<<" "<<bbox[1]<<" fsize dup add add sub\n";
        out<<
        "  SADDEL dup length 1 sub get 2 get /maxy exch def % max height\n"
        "  9999999 LEAF { aload pop exch pop min } forall\n"
        "  /miny exch def % min height\n"
        /* "  LEAF 0 get 1 get /miny exch def % min height\n" */
        "  maxy miny sub dup 20 div /MinHeight exch def\n"
        "  div scale\n"
        "  .5 LEAF 0 get 1 get neg translate\n"
        "  SetLabelFont\n"
        "  Drawlabels\n"
        "  DrawScale\n"
        "  SetBarFont\n"
        "  Connectlmins\n"
        "  showpage\n"
        "end\n"
        "%%EOF\n";

          // reset static pointer
    }

    bool
    LT_BarrierTree::cmp_saddle::operator()(int A, int B) 
    {
        assertbiu(leaves != NULL, "leaves variable not set");
        float diff = leaves->at(A).saddle_height - leaves->at(B).saddle_height;
        if (diff < -1e-6) 
            return true;
        else if (diff>1e-6) 
            return false;
        diff = leaves->at(A).height - leaves->at(B).height;
        return (diff<0)?true:false;
    }
      

}