main page namespaces classes files
ELL home

src/ell/rna/DW_RNA.cc

Go to the documentation of this file.

#include "ell/rna/DW_RNA.hh"

#include <biu/assertbiu.hh>
#include <biu/RandomNumberFactory.hh>

#include "ell/rna/S_RNAfe_SingleM.hh"
#include "ell/StateCollector.hh"

#include <algorithm>

namespace ell {

    DW_RNA_MorganHiggs::
	DW_RNA_MorganHiggs( const size_t iterations_)
     :  iterations(iterations_)
     {
     }
    
    DW_RNA_MorganHiggs::
	~DW_RNA_MorganHiggs()
    {
    }
    
    size_t
    DW_RNA_MorganHiggs::
	getIterations( void ) const
    {
        return iterations;
    }
    


    double 
    DW_RNA_MorganHiggs::
	addBond(   std::string& sequence
                , std::vector<size_t>& v_struct
                , size_t start
                , size_t end
                , StateCollector* sc) 
    
    {
        assertbiu(sc != NULL, "no StateCollector sc given (==NULL)");
        v_struct[start] = end;
        
        std::string structure(v_struct.size(),'.');
        
        for (size_t i = 0; i < v_struct.size(); ++i) {
            if(v_struct[i]!=0) {
                structure[i]='(';
                structure[v_struct[i]]=')';
            }
        }
        
        S_RNAfe_SingleM* rna = new S_RNAfe_SingleM(sequence,structure);
        double energy = rna->getEnergy();
        
        sc->add(*rna);
        
        delete rna;
        return energy;
    }

    double 
    DW_RNA_MorganHiggs::
	deleteBond(    std::string& sequence
                , std::vector<size_t>& v_struct
                , size_t start
                , StateCollector* sc) 
    
    {
        assertbiu(sc != NULL, "no StateCollector sc given (==NULL)");
        
        v_struct[start] = 0;
        std::string structure(v_struct.size(),'.');
            
        for (size_t i = 0; i < v_struct.size(); ++i) {
            if(v_struct[i]!=0) {
                structure[i]='(';
                structure[v_struct[i]]=')';
            }
        }
        
        S_RNAfe_SingleM* rna = new S_RNAfe_SingleM(sequence,structure);
        double energy = rna->getEnergy();
        
        sc->add(*rna);
        
        delete rna;
        return energy;
    }

    StateCollector*
    DW_RNA_MorganHiggs::
	walk(  const State* const begin
            , const State* const end
            , StateCollector* const scWalk
        ) const
    {
        using namespace std;
        
        typedef vector<size_t>::iterator vit;
        
    //  start = clock();
        assertbiu(begin != NULL, "State A is NULL");
        assertbiu(end != NULL, "State B is NULL");
        assertbiu(scWalk != NULL, "StateCollector is NULL");
        
        const S_RNAfe_SingleM* rna_A 
                    = dynamic_cast<const S_RNAfe_SingleM*>(begin);
        const S_RNAfe_SingleM* rna_B 
                    = dynamic_cast<const S_RNAfe_SingleM*>(end);
        
        assertbiu(rna_A != NULL, "state A is no RNA");
        assertbiu(rna_B != NULL, "state B is no RNA");

        assertbiu(rna_A->getSequence().compare(rna_B->getSequence()) == 0
                , "state A and B have different sequences");
        
        string sequence = rna_A->getSequence();
        string str_A = rna_A->getStructure();
        string str_B = rna_B->getStructure();
        size_t length = str_A.size();
        
        vector<size_t> structA, structB;
        vector<size_t> tempA, tempB;
            
        for (size_t i = 0; i < length; ++i) {
            structA.push_back(0);
            structB.push_back(0);
            
            if(str_A[i]=='(') tempA.push_back(i);
            if(str_A[i]==')') {
                structA[tempA.back()]=i;
                tempA.pop_back();
            }
                    
            if(str_B[i]=='(') tempB.push_back(i);
            if(str_B[i]==')') {
                structB[tempB.back()]=i;
                tempB.pop_back();
            }
        }
        
        vector<size_t> indexA, indexB;
        
        for (size_t i = 0; i < length; ++i) {
            if(structA[i] != structB[i]) {
                if(structA[i] != 0) indexA.push_back(i);
                if(structB[i] != 0) indexB.push_back(i);
            }
        }

        
        vector<vector<size_t> > conflicts;
        vector<size_t> noconflicts;
        vector<size_t> temp;
        vector<size_t> temp2(indexA.size(),0);
            
        for (size_t i = 0; i < indexB.size(); ++i) {
            temp.clear();
            for (size_t j = 0; j < indexA.size(); ++j) {
                if( !((indexB[i] < indexA[j] 
                       &&  structB[indexB[i]] > structA[indexA[j]]) 
                    || (indexB[i] > indexA[j] 
                       &&  structB[indexB[i]] < structA[indexA[j]]))) 
                {
                    temp.push_back(j);
                    temp2[j] = 1;
                }
            }       
            conflicts.push_back(temp);
        }
        
        for (size_t i = 0; i < indexA.size(); ++i) {
            if(temp2[i]==0) noconflicts.push_back(i); 
        }
        
        vector<size_t> choices;
        vector<size_t> addedBonds, deletedNonConflBonds;
        vector<size_t> temp_struct;
        vector<size_t> temp_struct_best;
        vector<vector<size_t> > temp_conflicts;
        vector<size_t> temp_noconflicts;
        size_t min_conflicts;
        size_t choice_add, choice_delete, toAdd, toDelete;
        SC_Listing *scl = NULL, *scl_best = NULL;
        double energy = 0, peakEnergy = -INT_MAX, Emax = INT_MAX;
        
    //  bucket1 = clock() - start;
        
        
        for (size_t count = 0; count < iterations; ++count) {
    //      cerr << "iteration "<<count<<endl;
            scl = new SC_Listing();
    //      scl->add(*rna_A);
            
            temp_struct = structA;
            temp_conflicts = conflicts;
            addedBonds.clear();
            peakEnergy = -INT_MAX;
    //      start = clock();
                    
            while(addedBonds.size() != conflicts.size()) {
                min_conflicts = INT_MAX;
                for (size_t i = 0; i < temp_conflicts.size(); ++i) {
                    if(find(addedBonds.begin(),addedBonds.end(),i) 
                            != addedBonds.end()) 
                    {
                        continue;
                    }
    //              if(count < 1) cerr << " bond "<<i<< " with size "<<temp_conflicts[i].size()<<endl;
                    if(temp_conflicts[i].size() < min_conflicts) { 
                        choices.clear(); 
                        min_conflicts = temp_conflicts[i].size();
                    } 
                    if(temp_conflicts[i].size() == min_conflicts) {
                        choices.push_back(i);
                    }
                }
    //          cerr << "choices to Add are:\t";
    //          for (size_t i = 0; i < choices.size(); ++i) {
    //              cerr << choices[i] << " ";              
    //          }
    //          cerr << endl;
                choice_add = (size_t)((double)(choices.size()) * 
                             (double)(biu::RNF::getRN(100)) / (double)100);
                toAdd = choices[choice_add];
    //          if(count<1) {cerr << "toAdd: " << toAdd<<endl;}
                
                while(temp_conflicts[toAdd].size() != 0) {
                    choice_delete = (size_t)(
                                    (double)(temp_conflicts[toAdd].size()) 
                                    * (double)(biu::RNF::getRN(100))
                                    / 100.0);
                    toDelete = temp_conflicts[toAdd][choice_delete];
                    energy = deleteBond(    sequence, temp_struct, 
                                            indexA[toDelete], scl);
                    if(energy > Emax) break;
                    if(energy > peakEnergy) peakEnergy = energy;
                    
    //              if(count<1) {cerr << "toDelete: " << toDelete<<endl;}
                    for (size_t i = 0; i < temp_conflicts.size(); ++i) {
                                            
                        for (size_t j = 0; j < temp_conflicts[i].size(); ++j) {
                            if(toDelete == temp_conflicts[i][j]) { 
                                temp_conflicts[i].erase(
                                        temp_conflicts[i].begin()+j);
                            }
                        }
    //                  if(count<1) {
    //                      cerr << "temp_conflicts["<<i<<"]:\t";
    //                      for (size_t k = 0; k < temp_conflicts[i].size(); ++k) {
    //                          cerr << temp_conflicts[i][k] << " ";
    //                      }
    //                      cerr << endl;
    //                  }
                    }
                }

                if(energy > Emax) break;
                            
                energy = addBond(   sequence, temp_struct, indexB[toAdd]
                                    , structB[indexB[toAdd]], scl);
                addedBonds.push_back(toAdd);
                
                if(energy > peakEnergy) peakEnergy = energy;
                if(energy > Emax) break;                
            }
    //      bucket2 += clock() - start;

            if(energy > Emax) { delete scl; continue;}
            
    //      start = clock();
            temp_noconflicts = noconflicts;
            deletedNonConflBonds.clear();
            while(temp_noconflicts.size() != 0) {
                choice_delete = (size_t)((double)(temp_noconflicts.size()) * 
                                (double)(biu::RNF::getRN(100)) / (double)100);
                toDelete = temp_noconflicts[choice_delete];
                energy = deleteBond(    sequence, temp_struct
                                        , indexA[toDelete], scl);
                if(energy > Emax) break;
                if(energy > peakEnergy) peakEnergy = energy;
                temp_noconflicts.erase( temp_noconflicts.begin() 
                                        + choice_delete);
            }
                
    //      bucket3 += clock() - start;
            if(energy > Emax) { delete scl; continue;}

            Emax = peakEnergy;
            delete scl_best;
            scl_best = scl;
            temp_struct_best = temp_struct;
        }
        
        list<State*> l = scl_best->getList();
        
        for(list<State*>::const_iterator cit = l.begin(); cit!=l.end(); ++cit)
        {
            scWalk->add(**cit);
        }
        

        string str_temp(temp_struct_best.size(),'.');
            
        for (size_t i = 0; i < temp_struct_best.size(); ++i) {
            if(temp_struct_best[i]!=0) {
                str_temp[i]='(';
                str_temp[temp_struct_best[i]]=')';
            }
        }
        
    //  cerr << "str_A:\t\t" << str_A<<endl;
    //  cerr << "str_temp:\t" << str_temp<<endl;
    //  cerr << "str_B:\t\t" << str_B<<endl;
        
        delete scl_best;
        
    //  cerr << "1: " << double(bucket1)/CLOCKS_PER_SEC<<endl;
    //  cerr << "2: " << double(bucket2)/CLOCKS_PER_SEC<<endl;
    //  cerr << "3: " << double(bucket3)/CLOCKS_PER_SEC<<endl;
        
        
        return scWalk;
    }


} // namespace ell