main page namespaces classes files
ELL home

src/ell/LT_SaddleNet.cc

Go to the documentation of this file.

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

#include <map>
#include <cmath>

namespace ell {

    const std::string LT_SaddleNet::LT_ID = "SN";
    
    
    
    double LT_SaddleNet::MISSING_SADDLE_PENALTY = 5.0;

    
    
    LT_SaddleNet::
	~LT_SaddleNet()
    {}


    // Calculates the relative 'distance' between landscape topologies.
    // On the LandscapeTopology level it is done on the minima only. For each
    // minimum m present in only one of the two LTs a value of 
    // MISSING_MINIMUM_PENALTY*(e^(-energy(m)/BOLTZMANN_KT)) is summed. The 
    // resulting sum is normalized by the number of (unique) minima in both
    // topologies.
    // For all O(n^2) saddles in the nets the RMSD is calculated and 
    // summed to the minima penalty. For each saddle present in one of the
    // saddle nets only, the highest saddle of both nets  
    // is used in the RMSD calculation for the missing saddles.
    // @param lt the saddle net to calculate the distance to
    // @return the distance
    double
    LT_SaddleNet::
	getDistance(const LandscapeTopology* const lt) const
    {
        assertbiu(lt!=NULL, "given LandscapeTopology to compare to is NULL!");

          // check if comparing with itself
        if (this == lt)
            return 0.0;
        
          // check if lt is an LT_MinimaSet object
        const LT_SaddleNet* const sn 
            = dynamic_cast<const LT_SaddleNet* const>(lt);
        if ( sn == NULL ) { 
            // no LT_SaddleNet --> use default function
            return LandscapeTopology::getDistance(lt);
        }

        double minimaDistance = 0.0;
        
          // sn1 will be the saddle net with the smaller number of minima
        const LT_SaddleNet* sn1 = this;
        const LT_SaddleNet* sn2 = sn;
        if (this->getMinCount() > sn->getMinCount()) {
            sn1 = sn;
            sn2 = this;
        }
            
        
        std::map<size_t, size_t> map2to1;


        size_t numOfDifferentMinima = 0;
        size_t curMatch = INVALID_INDEX;
          // find mapping of this minima to the minima of sn
        for (size_t i=sn1->getMinCount(); i-->0; ) {
              // try to map the current minimum of sn1
            curMatch = sn2->getMinIndex(sn1->getMin(i));
              // check if the minimum could be mapped
            if (curMatch != INVALID_INDEX) {
                  // store mapping
                map2to1[curMatch] = i;
            } else {
                  // add missing minimum penalty
                minimaDistance += MISSING_MINIMUM_PENALTY 
                                    * exp( -(sn1->getMin(i)->getEnergy())
                                            / BOLTZMANN_KT );
            }
        }
          // all minima of sn1 are considered to be unique
        numOfDifferentMinima += sn1->getMinCount();
          // find all minima of sn2 that could not be mapped
        for (size_t i=sn2->getMinCount(); i-->0; ) {
              // check if minimum of sn2 was not mapped
            if (map2to1.find(i) == map2to1.end()) {
                  // add missing minimum penalty
                minimaDistance += MISSING_MINIMUM_PENALTY 
                                    * exp( -(sn2->getMin(i)->getEnergy())
                                            / BOLTZMANN_KT );
                  // update number of different minima
                numOfDifferentMinima++;
            }
        }
          // normalize minima penalty by the number of different minima
        if (numOfDifferentMinima > 0) {
            minimaDistance /= (double)numOfDifferentMinima;
        }
        
        // saddle distance calculation
        
          // check all saddles
        double squareDeviation = 0.0;
        size_t numOfSaddles = 0;
        std::map<size_t,size_t>::const_iterator from = map2to1.begin();
        std::map<size_t,size_t>::const_iterator to = from;
        const State* saddleIn1 = NULL;
        const State* saddleIn2 = NULL;
        for ( from = map2to1.begin(); from != map2to1.end(); from++ ) {
            for ( ++(to = from); to != map2to1.end(); to++ ) {
                  // get the saddle states if existing
                saddleIn1 = sn1->getSaddle( from->second, to->second );
                saddleIn2 = sn2->getSaddle( from->first, to->first );
                  // check if there are saddles to compare 
                if ( saddleIn1 == NULL && saddleIn2 == NULL ) {
                    continue;
                } else if ( saddleIn1 != NULL && saddleIn2 != NULL ) {
                      // calculate square deviation
                    squareDeviation 
                        += (saddleIn1->getEnergy() - saddleIn2->getEnergy())
                            * (saddleIn1->getEnergy() - saddleIn2->getEnergy());
                } else if (saddleIn1 != NULL) {
                      // calculate penalty for missing saddle in sn2
                    squareDeviation += MISSING_SADDLE_PENALTY
                                        * exp( -(saddleIn1->getEnergy())
                                                / BOLTZMANN_KT );
                } else {
                      // calculate penalty for missing saddle in sn1
                    squareDeviation += MISSING_SADDLE_PENALTY
                                        * exp( -(saddleIn2->getEnergy())
                                                / BOLTZMANN_KT );
                }
                  // count the contribution to the saddle deviation
                numOfSaddles++;
            }
        }
        
          // check if a saddle was found
        if (numOfSaddles == 0) {
            return minimaDistance;
        }
          // return minima distance + RMSD on saddles (including penalty)
        return minimaDistance + sqrt( squareDeviation / (double)numOfSaddles );
    }

} // namespace ell


#include <sstream>
#include "ell/LT_MinimaSet.hh"

namespace ell {

    size_t LT_SaddleNet_AM::MatrixResizeAdd = 10; 
    const State* const LT_SaddleNet_AM::NULLpointer = NULL; 

    
    LT_SaddleNet_AM::
	~LT_SaddleNet_AM()
    {
        this->clear();
    }
    
    
    LT_SaddleNet_AM::
	LT_SaddleNet_AM( const size_t matrixDim )
      : adjacent(matrixDim, matrixDim, NULLpointer)
        , numOfMin(0)
        , minAccess()
        , mfeState(NULL)
        , sorted(true)
        , sortedEnd(0)
        , saddleESum(0.0)
        , saddleNum(0)
    {
        
    }
    
    
    LT_SaddleNet_AM::
	LT_SaddleNet_AM( const LT_SaddleNet_AM & toCopy )
      : adjacent(toCopy.numOfMin, toCopy.numOfMin, NULLpointer)
        , numOfMin(toCopy.numOfMin)
        , minAccess(toCopy.minAccess)
        , mfeState(NULL)
        , sorted(toCopy.sorted)
        , sortedEnd(toCopy.sortedEnd)
        , saddleESum(toCopy.saddleESum)
        , saddleNum(toCopy.saddleNum)
    {
          // copy all states and set mfeState
        for (size_t m1=0; m1<numOfMin; m1++) {
            for (size_t m2=m1; m2<numOfMin; m2++) {
                if (toCopy.adjacent[m1][m2] != NULLpointer) {
                      // make a copy of the stored state
                    this->adjacent[m1][m2] = toCopy.adjacent[m1][m2]->clone();
                      // create link if not minimum state on diagonal
                    if (m1 != m2) {
                        this->adjacent[m2][m1] = this->adjacent[m1][m2];
                    } else if (toCopy.adjacent[m1][m2] == toCopy.mfeState) {
                          // set mfe State
                        mfeState = this->adjacent[m1][m2];
                    }
                }
            }
        }
    }
    
    LT_SaddleNet_AM&
    LT_SaddleNet_AM::
    operator = (const LT_SaddleNet_AM & toCopy ) 
    
    {
        if (this != &toCopy) {
            this->numOfMin = toCopy.numOfMin;
            this->minAccess.resize( toCopy.minAccess.size(), 0 );
            std::copy(  toCopy.minAccess.begin()
                        , toCopy.minAccess.end()
                        , this->minAccess.begin() );
            this->sorted = toCopy.sorted;
            this->sortedEnd = toCopy.sortedEnd;
            this->saddleESum = toCopy.saddleESum;
            this->saddleNum = toCopy.saddleNum;
            
              // resize matrix if too small or too big
            if (    adjacent.numRows() < numOfMin 
                    || adjacent.numRows()-MatrixResizeAdd > numOfMin )
            {
                this->adjacent.resize(  this->numOfMin
                                        , this->numOfMin
                                        , NULLpointer);
            }
              // copy all states and set mfeState
            for (size_t m1=0; m1<numOfMin; m1++) {
                for (size_t m2=m1; m2<numOfMin; m2++) {
                    if (toCopy.adjacent[m1][m2] != NULLpointer) {
                          // make a copy of the stored state
                        this->adjacent[m1][m2] = toCopy.adjacent[m1][m2]->clone();
                          // create link if not minimum state on diagonal
                        if (m1 != m2) {
                            this->adjacent[m2][m1] = this->adjacent[m1][m2];
                        } else if (toCopy.adjacent[m1][m2] == toCopy.mfeState) {
                              // set mfe State
                            mfeState = this->adjacent[m1][m2];
                        }
                    } else {
                          // overwrite possible old entries
                        this->adjacent[m1][m2] = NULLpointer;
                        this->adjacent[m2][m1] = NULLpointer;
                    }
                }
            }
            
        }
        return *this;
    }
    
    
    void
    LT_SaddleNet_AM::
	clear()
    {
        assertbiu(numOfMin <= adjacent.numRows(), "numOfMin does not correlate with adjacency matrix (>size)");
        
          // delete all minima and saddles and reset pointer
        for (size_t from=0; from<numOfMin; from++) {
            for (size_t to=from; to<numOfMin; to++) {
                if (adjacent[from][to] != NULLpointer) {
                    delete adjacent[from][to];
                    adjacent[from][to] = NULLpointer;
                    adjacent[to][from] = NULLpointer;
                }
            }
        }
          // reset remaining information
        numOfMin = 0;
        minAccess.clear();
        mfeState = NULL;
        sorted = true;
        sortedEnd = 0;
        saddleESum = 0.0;
        saddleNum = 0;
    }
    

    const size_t 
    LT_SaddleNet_AM::
	addMin(const State* const s) 
    
    {
        assertbiu(s != NULL, "no minimum given (s is NULL)");
        assertbiu(getMinIndex(s) == LandscapeTopology::INVALID_INDEX
                , "given minimum is known and already part of the landscape");
        
        State* stateclone = s->clone();
        
          // check if matrix resize is neccessary
        if (numOfMin == adjacent.numRows()) {
              // resize matrix and initialize new elements with NULL 
            adjacent.resize(    numOfMin+MatrixResizeAdd
                                , numOfMin+MatrixResizeAdd
                                , NULLpointer);
        }
          // add minimum to adjacency matrix
        adjacent[numOfMin][numOfMin] = stateclone;
        minAccess.push_back(numOfMin);
          // increase minima counter
        numOfMin++;
        
          // check if new minimum is mfe
        if( mfeState == NULL || stateclone->operator<( *mfeState ) )
        {
            mfeState = stateclone;
        }
        
          // update sorted property
        sorted = false;

          // return getMinCount()-1 that is the index of the new minimum 
        return (numOfMin-1);
    }

    bool
    LT_SaddleNet_AM::
	addSaddle(const size_t m1, const size_t m2, const State* const s) 
    
    {
        assertbiu(m1<numOfMin, "m1 out of range");
        assertbiu(m2<numOfMin, "m2 out of range");
        assertbiu(m1!=m2, "m1 and m2 are equal --> saddle to itself");
        
        size_t from = m1, to = m2;
        if (m1>m2) {
            from = m2;
            to = m1;
        }
        
          // add saddle to adjacency list (but copy only once)
        if (adjacent[from][to] == NULLpointer) {
            adjacent[from][to] = s->clone();
            adjacent[to][from] = adjacent[from][to];
              // handle new saddle
            saddleNum++;
            saddleESum += s->getEnergy();
            return true;
        } else if (s->operator<( *adjacent[from][to] )) {
            // new saddle is lower than known one
              // update saddle sum by decrease by old saddle energy
            saddleESum -= adjacent[from][to]->getEnergy();
              // delete old saddle
            delete adjacent[from][to];
              // set new saddle in adjacency list (but copy only once)
            adjacent[from][to] = s->clone();
            adjacent[to][from] = adjacent[from][to];
              // update saddle sum by increase with new saddle energy
            saddleESum += s->getEnergy();
            return true;
        }
        
          // no saddle handling done and therefore no change
        return false;
    }

    bool
    LT_SaddleNet_AM::
	addSaddle(const State* const m1, const State* const m2, const State* const s) 
    
    {
        assertbiu( m1!=NULL, "no minimum m1 given (==NULL)");
        assertbiu( m2!=NULL, "no minimum m2 given (==NULL)");
        assertbiu( s!=NULL,  "no saddle s given (==NULL)");
        
          // get indices of the minima
        size_t im1 = getMinIndex( m1 );
        size_t im2 = getMinIndex( m2 );
        
        assertbiu( im1 != INVALID_INDEX, "m1 is no know minimum");
        assertbiu( im2 != INVALID_INDEX, "m2 is no know minimum");
        
          // perform adding via member function
        return addSaddle( im1, im2, s );
    }

    const State* const
    LT_SaddleNet_AM::
	getMin(size_t i) const 
    {
        assertbiu(i < numOfMin, "minimum index i out of bound!");
        
        return adjacent[minAccess[i]][minAccess[i]];
    }


    const size_t
    LT_SaddleNet_AM::
	getMinIndex(const State* const s) const
    {
        assertbiu(s != NULL, "given minimum State to get index of is NULL");

        for (size_t i = 0; i < numOfMin; ++i)
            if(adjacent[minAccess[i]][minAccess[i]]->operator==(*s)) 
                return i;
        
        return LandscapeTopology::INVALID_INDEX;
    }

    const size_t
    LT_SaddleNet_AM::
	getMinCount() const 
    {
        return numOfMin;
    }

    const State* const
    LT_SaddleNet_AM::
	getMFEState() const 
    {
        assertbiu(numOfMin != 0
                , "no minimum present and therefore no mfe state too");
        return mfeState;    
    }


    void
    LT_SaddleNet_AM::
	sort() 
    
    {
          // sort only if necessary
        if (sorted)
            return;
          // sort minima
        
        const State* curState = NULL;
        for (size_t curPos=1; curPos<numOfMin; curPos++) {
              // get minimum behind sorted range
            curState = getMin(minAccess[curPos]);
              // find position where to insert the current position
            size_t toInsert = curPos;
            while (toInsert>0) {
                if (getMin(minAccess[toInsert-1])->operator<(*curState)) {
                    break;
                }
                toInsert--;
            }
              // check if insertion is needed
            if (toInsert != curPos) {
                  // swap the two entries
                std::swap(minAccess[toInsert], minAccess[curPos]);
            }
        }
        
          // set sorted
        sorted = true;
    }


    
    bool
    LT_SaddleNet_AM::
	isSorted() const 
    {
        return sorted;
    }


    std::pair<int,std::string>
    LT_SaddleNet_AM::
	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 saddle net to add the saddles
        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_SaddleNet::LT_ID) != 0 
                    && lt_id.compare(LT_MinimaSet::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 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,std::string(""));
    }
    

    void
    LT_SaddleNet_AM::
	write( std::ostream& out, 
            const bool writeGraph ) const
    {

        const std::string dummyID = "ell::State";
        LandscapeTopology::writeHeader( out
                                        , LT_SaddleNet::LT_ID
                                        , numOfMin
                                        , (numOfMin==0?dummyID:adjacent[0][0]->getID()));
        
          // put all minima in with..
          // format : //MINIMUM index toString()
        for (size_t i = 0; i < numOfMin; ++i) {
            LandscapeTopology::writeMinimum( out, adjacent[i][i], i );
        }
        
          // the minima connected by the current saddle to print
        std::vector<size_t> minima;
          // the index of the next saddle to print
        size_t saddleIdx = numOfMin;
          // write all saddle
        for (size_t m1=0; m1<numOfMin; m1++) {
            minima.push_back(m1);
            for (size_t m2=m1+1; m2<numOfMin; m2++) {
                if (adjacent[m1][m2] != NULLpointer) {
                    minima.push_back(m2);
                      // write saddle to stream
                    LandscapeTopology::writeSaddle( out
                                                    , adjacent[m1][m2]
                                                    , saddleIdx
                                                    , minima );
                      // remind the printed saddle
                    saddleIdx++;
                    minima.pop_back();
                }
            }
            minima.pop_back();
        }
        
        if (writeGraph) {
            out << "\n" << "//GRAPHVIZ"<<"\n";
            
            out << "graph SADDLENET {"<<"\n";
            for (size_t i = 0; i < numOfMin; ++i) {
                out << "\t\tM"<<i<<";"<<"\n";           
            }
            size_t saddleIdx = numOfMin;
            for (size_t m1=0; m1<numOfMin; m1++) {
                for (size_t m2=m1+1; m2<numOfMin; m2++) {
                    if (adjacent[m1][m2] != NULLpointer) {
                          // write saddle to stream
                        out << "\t\tM"<<m1<<" -- M" <<m2 <<";"<<"\n";           
                          // remind the printed saddle
                        saddleIdx++;
                    }
                }
            }
            out << "label = \"\\n\\nLT_SaddleNet Diagram\\nwith "<<numOfMin<<" minima \";"<<"\n";
            out << "fontsize=20;"<<"\n";
            out << "} "<<"\n";
        }
        
          // flush the output stream
        out.flush();
        
    }
    
      /* Checks if two minima are connected via a direct saddle.
       * @param m1 the index of the first minimum
       * @param m2 the index of the second minimum
       * @return true if there is a direct saddle stored in the topology
       */
    bool
    LT_SaddleNet_AM::
	isSaddle( const size_t m1, const size_t m2 ) const
    {
        assertbiu(m1<numOfMin, "m1 out of range");
        assertbiu(m2<numOfMin, "m2 out of range");
    
        if (m1!=m2)
            return adjacent[m1][m2] != NULLpointer;
        
        return false;
    }
    
    
      /* Access to the direct saddle connecting two minima if existing.
       * @param m1 the index of the first minimum
       * @param m2 the index of the second minimum
       * @return the saddle State or NULL if no direct saddle exists
       */
    const State*
    LT_SaddleNet_AM::
	getSaddle( const size_t m1, const size_t m2 ) const
    {
        assertbiu(m1<numOfMin, "m1 out of range");
        assertbiu(m2<numOfMin, "m2 out of range");
        
        return adjacent[m1][m2];
    }
    
    double
    LT_SaddleNet_AM::
	getAvgSaddleEnergy( void ) const
    {
        if (saddleNum == 0)
            return 0.0;
          // calculate average saddle height
        return saddleESum / (double)saddleNum;
    }

    size_t
    LT_SaddleNet_AM::
	getSaddleNumber( void ) const
    {
        return saddleNum;
    }

    
} // namespace ell