main page namespaces classes files
ELL home

src/ell/LT_MinimaSet.cc

Go to the documentation of this file.

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

#include <biu/assertbiu.hh>

#include "ell/LT_MinimaSet.hh"


namespace ell {

const std::string LT_MinimaSet::LT_ID = "MS";

LT_MinimaSet::LT_MinimaSet()
: mfeIndex(LandscapeTopology::INVALID_INDEX), sorted(true)
{ }

LT_MinimaSet::LT_MinimaSet( const LT_MinimaSet& toCopy )
{
    this->operator =(toCopy);
}


LT_MinimaSet&
LT_MinimaSet::
operator = ( const LT_MinimaSet& toCopy )
{
      // check if toCopy is not the same object than this
    if (this == &toCopy) {
        return *this;
    }
      // clear data structures
    this->clear();
    
      // copy all minima
    vMinima.resize(toCopy.vMinima.size(), NULL);
    for (size_t i=0; i<vMinima.size(); i++) {
        vMinima[i] = toCopy.vMinima[i]->clone();
    }
      // copy remaining information
    mfeIndex = toCopy.mfeIndex;
    sorted = toCopy.sorted;

    return *this;
}

void 
LT_MinimaSet::clear() {
    for (size_t i = 0; i < vMinima.size(); ++i)
        delete vMinima[i];
    vMinima.clear();
    mfeIndex = LandscapeTopology::INVALID_INDEX;
    sorted = true;
}

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

bool
LT_MinimaSet::addSaddle(const size_t i, const size_t j, const State* const s) {
    // no saddle handling and therefore no change
    return false;
}

bool
LT_MinimaSet::addSaddle(const State* const si, const State* const sj, const State* const s) {
    // no saddle handling and therefore no change
    return false;
}

const size_t 
LT_MinimaSet::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");
    
    // add a copy of s to the set
    vMinima.push_back(s->clone());
    
     // 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 );
        }
    }

    return vMinima.size()-1;
}


const State* const
LT_MinimaSet::getMin(size_t i) const {
    assertbiu(i < vMinima.size(), "minimum index i out of bound!");
    
    return vMinima[i];
}


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

    for (size_t i = 0; i < vMinima.size(); ++i)
        if(*vMinima[i] == *s) 
            return i;
    
    return LandscapeTopology::INVALID_INDEX;
}

const size_t 
LT_MinimaSet::getMinCount() const {

    return vMinima.size();
}

const State* const 
LT_MinimaSet::getMFEState() const {
    assertbiu(mfeIndex != ell::LandscapeTopology::INVALID_INDEX, "no minimum present and therefore no mfe state too");
    return vMinima[mfeIndex];   
}

const std::vector<const State* > & 
LT_MinimaSet::getAllMin() const {
    return vMinima;
}



// 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 these not-matched minima.
// @param lt the landscape topology to calculate the distance to
// @return the distance
double
LT_MinimaSet::
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_MinimaSet* const ms = dynamic_cast<const LT_MinimaSet* const>(lt);
    if ( ms == NULL ) { 
        // no LT_MinimaSet --> use default function
        return LandscapeTopology::getDistance(lt);
    }

      // get access to all minima
    typedef std::vector<const State*> StateVec;
      // m1 are all sorted minima of this object
    StateVec* m1 = &(this->vMinima);
    if (!this->isSorted()) {
        m1 = new StateVec(vMinima);
        std::sort(m1->begin(), m1->end(), State::less);
    }
      // m2 are all sorted minima of ms object
    StateVec* m2 = &(ms->vMinima);
    if (!ms->isSorted()) {
        m2 = new StateVec(ms->vMinima);
        std::sort(m2->begin(), m2->end(), State::less);
    }
    
      // calculate distance via helper function
    double dist = LandscapeTopology::getSortedMinimaDistance(*m1,*m2);
    
      // garbage collection
    if ( ! this->isSorted() ) {
        delete m1;
    }
    if ( ! ms->isSorted() ) {
        delete m2;
    }
    
    return dist;
}


void 
LT_MinimaSet::sort() {
      // sort only if necessary
    if (sorted)
        return;
      // sort minima
    std::sort(vMinima.begin(),vMinima.end(),State::less);
      // update mfe index
    mfeIndex = 0;
      // set sorted
    sorted = true;
}

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

std::pair<int,std::string>
LT_MinimaSet::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 ) {
                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());
            }
        } 
        
          // 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_MinimaSet::write(    std::ostream& out, 
                        const bool writeGraph ) const
{
    size_t dim = vMinima.size();

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

      // put all minima in with..
      // format : //MINIMUM index toString()
    for (size_t i = 0; i < dim; ++i) {
        LandscapeTopology::writeMinimum( out, vMinima.at(i), i );
    }
    
    if (writeGraph) {
        out << "\n" << "//GRAPHVIZ"<<"\n";
        
        out << "graph SADDLENET {"<<"\n";
        for (size_t i = 0; i < dim; ++i) {
            out << "\t\tM"<<i<<";"<<"\n";           
        }   
        out << "label = \"\\n\\nLT_MinimaSet Diagram\\nwith "<<dim<<" minima \";"<<"\n";
        out << "fontsize=20;"<<"\n";
        out << "} "<<"\n";
    }
    
      // flush the output stream
    out.flush();
    
}


}//ell