main page namespaces classes files
ELL home

src/ell/LandscapeTopology.cc

Go to the documentation of this file.

#include "ell/LandscapeTopology.hh"

#include "biu/assertbiu.hh"

#include <sstream>
#include <vector>
#include <iomanip>
#include <cmath>
#include <algorithm>

namespace ell {


    const size_t 
    LandscapeTopology::
	INVALID_INDEX = UINT_MAX;
    
    // Avogadro constant N_A = 6.02214179 * 10^(23)
    // unit : 1/mol
    const double
    LandscapeTopology::
	AVOGADRO_CONSTANT_NA = 602214179000000000000000.0;
    
    // Gas constant R = 8.314472
    // unit : Joule / (Kelvin * mol)
    const double
    LandscapeTopology::
	GAS_CONSTANT_R = 8.314472;
    
    // Boltzmann constant  = R / N_A = 1.3806504 * 10^(-23) 
    // unit : Joule / Kelvin
    const double
    LandscapeTopology::
	BOLTZMANN_CONSTANT_KB = GAS_CONSTANT_R / AVOGADRO_CONSTANT_NA;
    
      // keyword used by the read and write method to identify information 
      // on the landscape topology type etc.
    const std::string 
    LandscapeTopology::
	OUTPUT_KEY_ELL_HEAD = std::string("ELL_LT");
    
      // keyword used by the read and write method to identify information 
      // on a minimum
    const std::string 
    LandscapeTopology::
	OUTPUT_KEY_MINIMUM = std::string("MINIMUM");
      // keyword used by the read and write method to identify information 
      // on a barrier
    const std::string 
    LandscapeTopology::
	OUTPUT_KEY_BARRIER = std::string("BARRIER");
      // keyword used by the read and write method to identify information 
      // on a saddle point
    const std::string 
    LandscapeTopology::
	OUTPUT_KEY_SADDLE = std::string("SADDLE");
    
    
    // temperature T = 37 Celsius = 310.15 Kelvin
    double 
    LandscapeTopology::
	BOLTZMANN_KT = BOLTZMANN_CONSTANT_KB * 310.15;

    double 
    LandscapeTopology::
	MISSING_MINIMUM_PENALTY = 1.0;
    

    LandscapeTopology::~LandscapeTopology() 
    
    {
    }

    
    // 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)) 
    // @param lt the landscape topology to calculate the distance to
    // @return the distance
    double
    LandscapeTopology::
	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;
        
          // get access to all minima
        typedef std::vector<const State*> StateVec;
        StateVec m1(this->getMinCount(),NULL);
        for (size_t i=0; i<m1.size(); i++) {
            m1[i] = this->getMin(i);
        }
        StateVec m2(lt->getMinCount(),NULL);
        for (size_t i=0; i<m2.size(); i++) {
            m2[i] = lt->getMin(i);
        }
          // ensure that ordered
        if (!this->isSorted()) {
            std::sort(m1.begin(), m1.end(), State::less);
        }
        if (!lt->isSorted()) {
            std::sort(m2.begin(), m2.end(), State::less);
        }
        
          // calculate distance via helper function
        return getSortedMinimaDistance(m1,m2);
    }
    
    double
    LandscapeTopology::
	getSortedMinimaDistance(   const std::vector<const State*> & m1,
                                const std::vector<const State*> & m2 )
    {
        double dist = 0.0;
        size_t numOfDifferentMinima = 0;
        
          // check if one set is empty
        if (m1.size() == 0 || m2.size() == 0) {
              // set m to the non-empty one
            const std::vector<const State*> & m = (m1.size() == 0?m2:m1);
              // check if both are empty
            if (m.size() == 0)
                return dist;
              // calc the penalty for all in the non-empty one
            for (size_t i=0; i<m.size(); i++) {
                dist += MISSING_MINIMUM_PENALTY 
                        * exp( -(m[i]->getEnergy()) / BOLTZMANN_KT );
            }
            return dist / (double)m.size();
        }
        
          // find all minima that are unique in one of the two topologies
        size_t i1 = 0;  // current index in m1 that is checked
        size_t i2 = 0;  // current index in m2 that is checked
        while (i1<m1.size() && i2<m2.size()) {
            numOfDifferentMinima++;
            if (*m1[i1] == *m2[i2]) {
                  // minimum in both topologies present --> go to next
                i1++;
                i2++;
            } else if ( State::less(m1[i1],m2[i2]) ) {
                  // smaller minimum only in m1 present --> handle distance
                dist += MISSING_MINIMUM_PENALTY 
                        * exp( -(m1[i1]->getEnergy()) / BOLTZMANN_KT );
                i1++;
            } else {
                  // smaller minimum only in m2 present --> handle distance
                dist += MISSING_MINIMUM_PENALTY 
                        * exp( -(m2[i2]->getEnergy()) / BOLTZMANN_KT );
                i2++;
            }
        }
          // handle remaining minima only present in m1
        while (i1 < m1.size()) {
            numOfDifferentMinima++;
            dist += MISSING_MINIMUM_PENALTY 
                    * exp( -(m1[i1]->getEnergy()) / BOLTZMANN_KT );
            i1++;
        }
          // handle remaining minima only present in m2
        while (i2 < m1.size()) {
            numOfDifferentMinima++;
            dist += MISSING_MINIMUM_PENALTY 
                    * exp( -(m2[i2]->getEnergy()) / BOLTZMANN_KT );
            i2++;
        }
          // normalize distance sum by number of different minima present in both
        if (numOfDifferentMinima == 0)
            return dist;
        else
            return dist / (double)numOfDifferentMinima;
    }
    
    void
    LandscapeTopology::
	writeHeader(   std::ostream & out, 
                    const std::string& LT_ID,
                    const size_t minNumber,
                    const std::string & StateDescription )
    {
        out <<"//" <<OUTPUT_KEY_ELL_HEAD <<"("
            <<LT_ID <<","
            <<minNumber <<","
            <<StateDescription 
            <<")\n";
    }
    
    std::pair<int,std::string>
    LandscapeTopology::
	parseHeader(   const std::string& line, 
                    std::string& LT_ID,
                    size_t& minNumber,
                    std::string & StateDescription )
    {
        typedef std::pair<int,std::string> RETURNTYPE;
        size_t start=0, end=0;
        std::stringstream sstr;
        
          // check if this line really encodes a barrier
        start = line.find(OUTPUT_KEY_ELL_HEAD);
        if(start == std::string::npos) {
            sstr.clear();
            sstr <<"no header description in given line";
            return RETURNTYPE(-3,sstr.str());
        }
          // jump over leading '(' of content
        start += 1+OUTPUT_KEY_ELL_HEAD.size();
            
          // read LT_ID index
        start = line.find_first_not_of(' ',start);
        end = line.find_first_of(',',start);
        if (start == end || start == std::string::npos) {
            sstr.clear();
            sstr <<OUTPUT_KEY_ELL_HEAD <<" found but no LT_ID given";
            return RETURNTYPE(-3,sstr.str());
        }
        LT_ID = line.substr(start, end-start);
        
          // read minNumber
        start = line.find_first_not_of(' ',end+1);
        end = line.find_first_of(',',start);
        if (start == end || start == std::string::npos) {
            sstr.clear();
            sstr <<OUTPUT_KEY_ELL_HEAD
                <<" found but no number of minima given";
            return RETURNTYPE(-1,sstr.str());
        }
        size_t num = INT_MAX;
        sstr.clear(); 
        sstr.str(line.substr(start, end-start));
        sstr >> num;
        if (num == INT_MAX) {
            sstr.clear();
            sstr <<OUTPUT_KEY_ELL_HEAD 
                <<" found but minima number not parsable into integer";
            return RETURNTYPE(-1,sstr.str());
        }
        if (num < 0) {
            sstr.clear();
            sstr <<OUTPUT_KEY_ELL_HEAD 
                <<" found but minima number smaller than zero";
            return RETURNTYPE(-1,sstr.str());
        }
        minNumber = num;
        
          // read StateDescriptor
        start = line.find_first_not_of(' ',end+1);
        end = line.find_last_of(')');
        end = line.find_last_not_of(' ', end);
        if (start >= end) {
            sstr.clear();
            sstr <<OUTPUT_KEY_ELL_HEAD 
                <<" found but no state description given";
            return RETURNTYPE(-1,sstr.str());
        }
        StateDescription = line.substr(start, end-start);

          // no error occured so far --> return new connection state s 
        return RETURNTYPE(0,std::string(""));
    }
    
    void
    LandscapeTopology::
	writeMinimum(  std::ostream & out, 
                    const State* const minimum,
                    const size_t index )
    {
        out << "//" <<OUTPUT_KEY_MINIMUM <<' '
            <<std::setw(7) <<index 
            <<std::setw(3)<<' '
            <<minimum->toString()
            <<"\n";
    }

    
    void
    LandscapeTopology::
	writeMinimum(  std::ostream & out, 
                    const std::vector<const State*>& minima,
                    const size_t index )
    {
        out << "//" <<OUTPUT_KEY_MINIMUM <<' '
            <<std::setw(7) <<index 
            <<std::setw(3)<<' ';
         // write all minima to stream
        for (size_t i=0; i<minima.size(); i++) {
            out <<minima[i]->toString() <<' ';
        }
        out <<"\n";
    }


    std::pair<int,std::string>
    LandscapeTopology::
	parseMinimum(  LandscapeTopology & toFill,
                    const std::string & line, 
                    const State * const templateState,
                    Int2SizeT_MAP & idx2min  )
    {
        typedef std::pair<int,std::string> RETURNTYPE;
        size_t start=0, end=0;
        State* s = NULL;
        std::stringstream sstr;
        
          // check if this line really encodes a barrier
        start = line.find(OUTPUT_KEY_MINIMUM);
        if(start == std::string::npos) {
            sstr.clear();
            sstr <<"no minimum description";
            return RETURNTYPE(-3,sstr.str());
        }
            
          // read minimum index
        start = line.find_first_not_of(' ',start+7);
        if (start == std::string::npos) {
            sstr.clear();
            sstr <<OUTPUT_KEY_MINIMUM <<" found but no index given";
            return RETURNTYPE(-3,sstr.str());
        }
        end = line.find_first_of(' ',start);
        int idx = INT_MAX;
        sstr.clear(); 
        sstr.str(line.substr(start, end-start));
        sstr >> idx;
        if (idx == INT_MAX) {
            sstr.clear();
            sstr <<OUTPUT_KEY_MINIMUM 
                <<" found but index not parsable into integer";
            return RETURNTYPE(-3,sstr.str());
        }
        if (idx2min.find(idx) != idx2min.end()) {
            sstr.clear();
            sstr <<OUTPUT_KEY_MINIMUM <<" index '"
                <<idx
                <<"' was used before (not unique)";
            return RETURNTYPE(-3,sstr.str());
        }
        
          // read first string representation from minima state list
          // and add to tree
        start = line.find_first_not_of(' ',end);
        if (start == std::string::npos) {
            sstr.clear();
            sstr <<OUTPUT_KEY_MINIMUM 
                <<" found but no state string representation given";
            return RETURNTYPE(-3,sstr.str());
        }
        end = line.find_first_of(' ',start+1);
        
          // create new state
        s = templateState->fromString( line.substr(start, end-start) );
        if(s == NULL) {
                sstr.clear();
                sstr <<OUTPUT_KEY_MINIMUM 
                    <<" state could not be generated from string representation";
            return RETURNTYPE(-3,sstr.str());
        }
          // add the new minimum to the landscape and store the mapping
        idx2min[idx] = toFill.addMin(s);
        
          // garbage collection
        if (s != NULL) {
            delete s;
            s = NULL;
        }
          // no error occured so far, return remaining line content
        if (end < line.size()) {
            return RETURNTYPE(0,line.substr(end));
        } else {
            return RETURNTYPE(0,std::string(""));
        }
    }
    
    void
    LandscapeTopology::
	writeBarrier(  std::ostream & out, 
                    const State* const barrier,
                    const size_t index,
                    const std::vector<size_t> & minConnected )
    {
        writeConnection(    out
                            , barrier
                            , index
                            , minConnected
                            , OUTPUT_KEY_BARRIER );
    }

    void
    LandscapeTopology::
	writeSaddle(   std::ostream & out, 
                    const State* const saddle,
                    const size_t index,
                    const std::vector<size_t> & minConnected )
    {
        writeConnection(    out
                            , saddle
                            , index
                            , minConnected
                            , OUTPUT_KEY_SADDLE );
    }

    void
    LandscapeTopology::
	writeConnection(   std::ostream & out, 
                        const State* const conState,
                        const size_t index,
                        const std::vector<size_t> & minConnected,
                        const std::string& connectionKey )
    {
        assertbiu(conState!=NULL, "no connection state given (==NULL)");
        
          // print connection index
        out << "//" <<connectionKey<<" " <<std::setw(7) <<index
            <<std::setw(3)<<" ";
          // print connected minima
        for (std::vector<size_t>::const_iterator it = minConnected.begin();
                it != minConnected.end(); it++) 
        {
            out <<(*it) <<" ";
        }
          // print connection object
        out <<conState->toString()
            <<"\n";

    }


    std::pair<int,std::string>
    LandscapeTopology::
	parseBarrier(  LandscapeTopology & toFill,
                    const std::string & line, 
                    const State * const templateState,
                    const Int2SizeT_MAP & idx2min  )
    {
        std::vector<size_t> minima;
        
          // forward parsing to non-specific member
        std::pair<State*,std::string> parseResult
            = parseConnection(  line
                                , templateState
                                , idx2min 
                                , OUTPUT_KEY_BARRIER
                                , minima );
        
          // check if a parsing error occured
        if (parseResult.first == NULL) {
            return std::pair<int,std::string>(-3, parseResult.second);
        }
        
          // add barrier to the landscape (only FIRST to ALL)
        size_t m0 = idx2min.find(minima[0])->second; 
        for (size_t m=1; m<minima.size(); m++) {
            toFill.addSaddle(   m0
                                , idx2min.find(minima[m])->second
                                , parseResult.first );
        }

          // parsing was successfull
        return std::pair<int,std::string>(0,std::string(""));
    }
        
    std::pair<int,std::string>
    LandscapeTopology::
	parseSaddle(   LandscapeTopology & toFill,
                    const std::string & line, 
                    const State * const templateState,
                    const Int2SizeT_MAP & idx2min  )
    {
        
        std::vector<size_t> minima;
        
          // forward parsing to non-specific member
        std::pair<State*,std::string> parseResult 
            = parseConnection(  line
                                , templateState
                                , idx2min 
                                , OUTPUT_KEY_SADDLE
                                , minima );
        
          // check if a parsing error occured
        if (parseResult.first == NULL) {
            return std::pair<int,std::string>(-3, parseResult.second);
        }
        
          // add the saddle to the landscape (all to all)
        for (size_t m1=0; m1<minima.size(); m1++) {
            for (size_t m2=m1+1; m2<minima.size(); m2++) {
                toFill.addSaddle(   idx2min.find(minima[m1])->second
                                    , idx2min.find(minima[m2])->second
                                    , parseResult.first );
            }
        }

          // parsing was successfull
        return std::pair<int,std::string>(0,std::string(""));
    }
        

    std::pair<State*,std::string>
    LandscapeTopology::
	parseConnection(   const std::string & line, 
                        const State * const templateState,
                        const Int2SizeT_MAP & idx2min,
                        const std::string& connectionKey,
                        std::vector<size_t> & minima )
    {
        
        typedef std::pair<State*,std::string> RETURNTYPE;
        size_t start=0, end=0;
        State* s = NULL;
        std::stringstream sstr;
        
          // clear output parameter where the minima are stored 
        minima.clear();
        
          // check if this line really encodes a connection of the given type
        start = line.find(connectionKey);
        if(start == std::string::npos) {
            sstr.clear();
            sstr <<"no " <<connectionKey
                <<" connection description present";
            return RETURNTYPE(NULL,sstr.str());
        }
        
          // read connection index
        start = line.find_first_not_of(' ',start+connectionKey.size());
        if (start == std::string::npos) {
            sstr.clear();
            sstr <<connectionKey
                <<" found but no connection index given";
            return RETURNTYPE(NULL,sstr.str());
        }
        end = line.find_first_of(' ',start);
        size_t idx = INT_MAX;
        sstr.clear(); 
        sstr.str(line.substr(start, end-start));
        sstr >> idx;
        if (idx == INT_MAX) {
            sstr.clear();
            sstr <<connectionKey 
                <<" found but connection index not parsable into integer";
            return RETURNTYPE(NULL,sstr.str());
        }
        
          // find substring that contains connection State string representation
        size_t s_end = line.find_last_not_of(' ')+1;
        if (s_end < end) {
            sstr.clear();
            sstr <<connectionKey 
                <<" found but no state string representation given";
            return RETURNTYPE(NULL,sstr.str());
        }
        size_t s_start = line.find_last_of(' ', s_end -1) +1;
        
          // read first string representation and add to tree
          // create new state
        s = templateState->fromString( line.substr(s_start, s_end-s_start) );
        if(s == NULL) {
                sstr.clear();
                sstr <<connectionKey 
                    <<" state could not be generated from string representation";
            return RETURNTYPE(NULL,sstr.str());
        }
        
          // get all minima to add barriers for
          // get string limits that hold all minima
        start = line.find_first_not_of(' ', end);
        end = line.find_last_not_of(' ', s_start -1)+1;
        if (start >= end) {
                sstr.clear();
                sstr <<connectionKey 
                    <<" found but no minima that are connected given";
                delete s;
                return RETURNTYPE(NULL,sstr.str());
        }
        sstr.clear();
        sstr.str(line.substr(start, end-start));
          // read all minima
        while (!sstr.fail()) {
            idx = INT_MAX;
            sstr >>idx;
            if (idx != INT_MAX) {
                if (idx2min.find(idx) == idx2min.end()) {
                        sstr.clear();
                        sstr <<connectionKey 
                            <<" found but connected minimum '"
                            <<idx <<"' not known so far";
                        delete s;
                        return RETURNTYPE(NULL,sstr.str());
                }
                minima.push_back(idx);
            } else {
                break;
            }
        }
        if (minima.size() == 0) {
                sstr.clear();
                sstr <<connectionKey 
                    <<" found but minima indices are not parsable";
                delete s;
                return RETURNTYPE(NULL,sstr.str());
        } else if (minima.size() == 1) {
                sstr.clear();
                sstr <<connectionKey 
                    <<" found but only ONE minimum index was parsable";
                delete s;
                return RETURNTYPE(NULL,sstr.str());
        }
        
          // no error occured so far --> return new connection state s 
        return RETURNTYPE(s,std::string(""));
    }

} // namespace ell