main page namespaces classes files
ELL home

src/ell/LT_MinimaSet.cc

Go to the documentation of this file.

#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)
{ }

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 (LandscapeTopology::isSmaller_E_S( 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 = LandscapeTopology::isSmaller_E_S( 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 RMSD between this and a given LT_MinimaSet(difference
double 
LT_MinimaSet::getRMSD(const LandscapeTopology* const pLT ) const {
    assertbiu(pLT!=NULL, "given LandscapeTopology is NULL!");
    
    const LT_MinimaSet* const ms = dynamic_cast<const LT_MinimaSet* const>(pLT);
    assertbiu(ms!=NULL, "given LandscapeTopology no LT_MinimaSet!");
        
    double rmsd = 0;
    
    
    std::vector<const State*> vM1 = vMinima;
    std::vector<const State*> vM2 = ms->vMinima;
    
    std::sort(vM1.begin(), vM1.end(), LandscapeTopology::isSmaller_E_S);
    std::sort(vM2.begin(), vM2.end(), LandscapeTopology::isSmaller_E_S);
    
    size_t i = 0;
    for (; i < vM1.size() && i < vM2.size(); ++i)
        rmsd += pow( vM1[i]->getEnergy() - vM2[i]->getEnergy(), 2);

    for (; i < vM1.size() || i < vM2.size(); ++i)
        if(i < vM1.size())
            rmsd += pow( vM1[i]->getEnergy() - vM2[0]->getEnergy(), 2);
        else // (i < vM2.size())
            rmsd += pow( vM1[0]->getEnergy() - vM2[i]->getEnergy(), 2);
    
    rmsd /= static_cast<double>(std::max(vM1.size(),vM2.size()));
    rmsd = sqrt(rmsd);
    
    return rmsd;
}

void 
LT_MinimaSet::sort() {
      // sort only if necessary
    if (sorted)
        return;
      // sort minima
    std::sort(vMinima.begin(),vMinima.end(),LandscapeTopology::isSmaller_E_S);
      // 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, 
                    const std::string & StateDescription )
{
    typedef std::pair<int,std::string> RETURNTYPE;
    State* s;
    enum STATUS {NOTHING, MINIMUM};
    STATUS status;
    
    std::string sequence, structure, line, substring;
    std::stringstream sstr;
    
    size_t type= INT_MAX, dim = INT_MAX;
    size_t index_i, index_j;
    size_t lineNumber = 0;
    
    std::string::size_type i, k, l;
    
    
    std::getline(in, line);
    
    i = line.find("ELL-LT-TYPE=");
    if(i != std::string::npos) { 
        line=line.substr(i+12);
        
        if (line.compare(LT_MinimaSet::LT_ID) != 0) {
            return RETURNTYPE(-1,"ELL-LT-TYPE not supported");
        }
    } else  {
        return RETURNTYPE(-1,"ELL-LT-TYPE missing in first line of file");
    }
    
    std::getline(in, line);
    
    i = line.find("STATETYPE=");
    if(i != std::string::npos) { 
        line=line.substr(i+10);
        
        if(line.compare(StateDescription) != 0) {
            RETURNTYPE(-2,"STATETYPE differs given type");
        } 
    } else  {
        RETURNTYPE(-2,"STATETYPE missing in second line of file");
    }

    std::getline(in, line);
    i = line.find("DIMENSION=");
    if(i != std::string::npos) {
        line=line.substr(i+10);         
        
        sstr.clear();
        sstr.str(line);
        sstr >> dim;
    } else  {
        RETURNTYPE(-3,"DIMENSION missing in third line of file");
    }
    
    
    lineNumber = 3;
    while( std::getline(in, line) ) {
        ++lineNumber;
        s = NULL;
        
//      cout << "line "<<lineNumber<<" is "<<line <<endl;
        index_i = index_j = INT_MAX;
        i = k = l = INT_MAX;
        
        status = NOTHING;
        
        i = line.find("GRAPHVIZ");
        if(i != std::string::npos) break;
        
        i = line.find("MINIMUM");
        if(i != std::string::npos) {
            if(type== INT_MAX || dim==INT_MAX) {
                RETURNTYPE(-3,"Wrong file format: Statetype or Dimension not defined before first state definition !");
            }
              // read minimum index
            size_t pos_i_begin = line.find("i=")+2;
            size_t pos_i_end   = line.find_first_of(" \t",pos_i_begin);
            sstr.clear();
            sstr.str( line.substr( pos_i_begin, pos_i_end ) );
            sstr >> index_i;
            
            status = MINIMUM;
            
            s = templateState->fromString(line.substr(line.find_first_not_of(" \t",pos_i_end)));
            if(s == NULL) {
                std::stringstream errstr;
                errstr <<"State object could not be generated from line "<<lineNumber <<" !";
                return RETURNTYPE(-3,errstr.str());
                continue;
            }
        } 
        
        if(status==NOTHING) continue;
 
        if(status==MINIMUM) 
            addMin(s);  // s is a new object already
        
        delete s;
        
    }

    
    return RETURNTYPE(0,"");
}

void
LT_MinimaSet::write(    std::ostream& out, 
                    const std::string & StateDescription ) const
{
    size_t dim = vMinima.size();


    out << "//ELL-LT-TYPE=" <<LT_ID <<"\n";
    out << "//STATETYPE=" <<StateDescription<< "\n";
    out << "//DIMENSION="<<dim<< "\n";
    
    for (size_t i = 0; i < dim; ++i) {
        out << "//MINIMUM\ti="<<i<<"\t"<<vMinima.at(i)->toString()<<"\n";
    }
    
    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