src/ell/rna/S_RNAfe_SingleM.cc
Go to the documentation of this file.00001 #include "ell/rna/S_RNAfe_SingleM.hh"
00002 #include <cmath>
00003 #include <biu/assertbiu.hh>
00004 #include "biu/RandomNumberFactory.hh"
00005
00006
00007 namespace ell
00008 {
00009
00010
00011 const std::string S_RNAfe_SingleM::ID = std::string("ell::S_RNAfe_SingleM");
00012
00013
00014 const size_t
00015 S_RNAfe_SingleM::RandomNeighborList::ItState::SWITCH_MODE_THRESHOLD = 20;
00016
00018
00020
00021 S_RNAfe_SingleM::S_RNAfe_SingleM(
00022 const std::string& rnaSeqStr,
00023 const std::string& rnaStructBracketNotStr)
00024 : rnaFreeEnergy(rnaSeqStr, rnaStructBracketNotStr)
00025 {
00026 assertbiu(rnaFreeEnergy.isValid(),
00027 "S_RNAfe_SingleM has to be a valid state.");
00028
00029 }
00030
00031 S_RNAfe_SingleM::S_RNAfe_SingleM(const RNAFreeEnergy& rna)
00032 : rnaFreeEnergy(rna)
00033 {
00034 assertbiu(rnaFreeEnergy.isValid(),
00035 "S_RNAfe_SingleM has to be a valid state.");
00036
00037 }
00038
00039 S_RNAfe_SingleM::S_RNAfe_SingleM(
00040 const S_RNAfe_SingleM& rnaFESMS)
00041 : rnaFreeEnergy(rnaFESMS.rnaFreeEnergy)
00042 {
00043 }
00044
00045 S_RNAfe_SingleM::~S_RNAfe_SingleM() {
00046 }
00047
00048
00049
00050
00051
00052 const std::string&
00053 S_RNAfe_SingleM::getID( void ) const {
00054 return S_RNAfe_SingleM::ID;
00055 }
00056
00057
00058 S_RNAfe_SingleM::SingleMove
00059 S_RNAfe_SingleM::firstValidSingleMove() const {
00060 return searchForNextValidSingleMove(0,1);
00061 }
00062
00063 S_RNAfe_SingleM::SingleMove
00064 S_RNAfe_SingleM::nextValidSingleMove(
00065 const SingleMove lastSingleMove) const
00066 {
00067
00068
00069
00070 if (rnaFreeEnergy.getStructureRef().at(lastSingleMove.first) != biu::RNAStructure::STRUCT_UNBOUND ) {
00071 return searchForNextValidSingleMove(lastSingleMove.first+1,
00072 lastSingleMove.first+2);
00073 }
00074 else {
00075 return searchForNextValidSingleMove(lastSingleMove.first,
00076 lastSingleMove.second+1);
00077 }
00078 }
00079
00080 S_RNAfe_SingleM::SingleMove
00081 S_RNAfe_SingleM:: searchForNextValidSingleMove( size_t pos1,
00082 size_t pos2) const
00083 {
00084 assertbiu(pos1 < pos2,
00085 "RNA single move search pos 1 has to be less than pos 2");
00086
00087
00088 if (pos2==rnaFreeEnergy.getLength() && pos1 == 0) {
00089 pos1++;
00090 pos2 = pos1+1;
00091 }
00092
00093
00094
00095
00096
00097
00098 while ( (pos1+rnaFreeEnergy.getMinLoopLength())
00099 < rnaFreeEnergy.getLength())
00100 {
00101 if ( rnaFreeEnergy.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_UNBOUND )
00102 {
00103 for(pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2-1);
00104 pos2 < rnaFreeEnergy.getLength();
00105 pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2))
00106 {
00107 if ((pos2 - pos1) > rnaFreeEnergy.getMinLoopLength() &&
00108 rnaFreeEnergy.isAllowedBasePair(pos1, pos2)) {
00109
00110
00111 return SingleMove(pos1,pos2);
00112 break;
00113 }
00114 }
00115 }
00116 else if ( rnaFreeEnergy.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_BND_OPEN )
00117 {
00118 pos2 = rnaFreeEnergy.getClosingBond(pos1);
00119
00120 return SingleMove(pos1,pos2);
00121 }
00122
00123 pos1++;
00124 pos2 = pos1;
00125 }
00126
00127 return SingleMove(0,biu::RNAStructure::INVALID_INDEX);
00128 }
00129
00130 bool
00131 S_RNAfe_SingleM::operator== (const State& state2) const {
00132 return rnaFreeEnergy == (dynamic_cast<
00133 const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
00134 }
00135
00136 bool
00137 S_RNAfe_SingleM::operator!= (const State& state2) const {
00138 return rnaFreeEnergy != (dynamic_cast<
00139 const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
00140 }
00141
00142 bool
00143 S_RNAfe_SingleM::operator< (const State& rna2) const
00144 {
00145 if (this == &rna2) {
00146 return false;
00147 }
00148 const S_RNAfe_SingleM* s2 = dynamic_cast<const S_RNAfe_SingleM*>(&rna2);
00149 assertbiu(s2!=NULL, "rna2 is no S_RNAfe_SingleM object");
00150 assertbiu(rnaFreeEnergy.getLength() == s2->rnaFreeEnergy.getLength()
00151 , "the two RNA molecules differ in length");
00152
00153 return (this->getEnergy() < s2->getEnergy())
00154 || (this->getEnergy() == s2->getEnergy()
00155 && std::lexicographical_compare(
00156 this->rnaFreeEnergy.getStructureRef().begin()
00157 , this->rnaFreeEnergy.getStructureRef().end()
00158 , s2->rnaFreeEnergy.getStructureRef().begin()
00159 , s2->rnaFreeEnergy.getStructureRef().end()) );
00160 }
00161
00162 double
00163 S_RNAfe_SingleM::getEnergy(void) const {
00164 return rnaFreeEnergy.getEnergy();
00165 }
00166
00167 unsigned int
00168 S_RNAfe_SingleM::getMinimalDistance(
00169 const State& _state2) const {
00170 S_RNAfe_SingleM state2 =
00171 dynamic_cast<const S_RNAfe_SingleM&>(_state2);
00172
00173 const biu::Structure& state1Struct = rnaFreeEnergy.getStructureRef();
00174 const biu::Structure& state2Struct = state2.rnaFreeEnergy.getStructureRef();
00175
00176 assertbiu(state2Struct.size() == state1Struct.size(),
00177 "Length of both RNA has to be equal to compute the distance.");
00178
00179 unsigned int distance = 0;
00180 for (size_t i = 0; i < state1Struct.size(); i++){
00181
00182 if (state1Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN &&
00183 state2Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN ) {
00184
00185 if (rnaFreeEnergy.getClosingBond(i) !=
00186 state2.rnaFreeEnergy.getClosingBond(i)) {
00187 distance += 2;
00188 }
00189 }
00190
00191 else if (state1Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN ||
00192 state2Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN ) {
00193 distance += 1;
00194 }
00195 }
00196 return distance;
00197 }
00198
00199 S_RNAfe_SingleM*
00200 S_RNAfe_SingleM::clone(State* toFill) const {
00201 if (toFill == NULL) {
00202 return new S_RNAfe_SingleM(*this);
00203 }
00204
00205 assertbiu( dynamic_cast<S_RNAfe_SingleM*>(toFill) != NULL
00206 , "given State is no S_RNAfe_SingleM instance");
00207 S_RNAfe_SingleM* s = static_cast<S_RNAfe_SingleM*>(toFill);
00208
00209
00210 s->rnaFreeEnergy = this->rnaFreeEnergy;
00211
00212 return s;
00213 }
00214
00215 State*
00216 S_RNAfe_SingleM::fromString(const std::string& stringRep) const {
00217
00218 std::string backboneRep = ".()";
00219 size_t pos = stringRep.find_first_not_of(backboneRep);
00220 const std::string rnaStructBracketDotStr = stringRep.substr(0, pos-1);
00221 assertbiu(rnaStructBracketDotStr.size() != 0, "string does not contain a dot-bracket structure encoding");
00222 assertbiu(rnaStructBracketDotStr.size() == this->rnaFreeEnergy.getLength(), "string is too short for this RNA molecule");
00223
00224 S_RNAfe_SingleM* ret = this->clone();
00225
00226 ret->rnaFreeEnergy.setStructure(
00227 ret->rnaFreeEnergy.getStructureAlphabet()->getSequence(
00228 rnaStructBracketDotStr));
00229 return ret;
00230 }
00231
00232 std::string
00233 S_RNAfe_SingleM::toString() const{
00234 return rnaFreeEnergy.getStructureString();
00235 }
00236
00237 std::string&
00238 S_RNAfe_SingleM::toString( std::string & toFill ) const {
00239
00240 toFill = rnaFreeEnergy.getStructureString();
00241 return toFill;
00242 }
00243
00244
00245 std::string
00246 S_RNAfe_SingleM::getSequence() const {
00247 return rnaFreeEnergy.getSequenceString();
00248 }
00249
00250 std::string
00251 S_RNAfe_SingleM::getStructure() const {
00252 return rnaFreeEnergy.getStructureString();
00253 }
00254
00255 State::NeighborListPtr
00256 S_RNAfe_SingleM::getNeighborList() const {
00257 return NeighborListPtr(new NeighborList(this));
00258 }
00259
00260 State::NeighborListPtr
00261 S_RNAfe_SingleM::getRandomNeighborList() const
00262 {
00263 return NeighborListPtr(new RandomNeighborList(this));
00264 }
00265
00266
00267 S_RNAfe_SingleM::NeighborList::NeighborList(
00268 const S_RNAfe_SingleM* _origin)
00269 : origin(_origin)
00270 {
00271 }
00272
00273 S_RNAfe_SingleM::NeighborList::~NeighborList() {
00274 }
00275
00276 State*
00277 S_RNAfe_SingleM::NeighborList::first(
00278 State::NeighborList::ItState** itstate) const {
00279 if (*itstate != NULL) delete *itstate;
00280
00281 ItState* myitstate = new ItState();
00282 *itstate = myitstate;
00283
00284 S_RNAfe_SingleM* rnaState = NULL;
00285
00286 SingleMove firstMove = origin->firstValidSingleMove();
00287
00288
00289 if (firstMove != SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00290 rnaState = new S_RNAfe_SingleM(*origin);
00291 rnaState->rnaFreeEnergy.applySingleMoveInPlace(firstMove);
00292 myitstate->lastMove = firstMove;
00293 }
00294 return rnaState;
00295 }
00296
00297 State*
00298 S_RNAfe_SingleM::NeighborList::next(
00299 State::NeighborList::ItState* itstate, State* elem) const {
00300 assertbiu(elem != NULL, "Elem is not allowed to be NULL");
00301
00302 ItState* myitstate = static_cast<ItState*>(itstate);
00303
00304
00305 SingleMove nextMove = origin->nextValidSingleMove(myitstate->lastMove);
00306
00307 if (nextMove == SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00308 return NULL;
00309 }
00310 else {
00311 S_RNAfe_SingleM* rnaState =
00312 dynamic_cast<S_RNAfe_SingleM*>(elem);
00313
00314 rnaState->rnaFreeEnergy.applySingleMoveInPlace(myitstate->lastMove);
00315
00316 rnaState->rnaFreeEnergy.applySingleMoveInPlace(nextMove);
00317 myitstate->lastMove = nextMove;
00318 return rnaState;
00319 }
00320 }
00321
00322
00323 S_RNAfe_SingleM::RandomNeighborList::RandomNeighborList(
00324 const S_RNAfe_SingleM* _origin)
00325 : origin(_origin)
00326 {
00327 }
00328
00329 S_RNAfe_SingleM::RandomNeighborList::~RandomNeighborList() {
00330 }
00331
00332 void
00333 S_RNAfe_SingleM::RandomNeighborList::switchMode(
00334 ItState* itstate) const
00335 {
00336
00337 ItState* myitstate = static_cast<ItState*>(itstate);
00338
00339 SingleMove nextMove = origin->firstValidSingleMove();
00340
00341
00342 while (nextMove != SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00343
00344 if (myitstate->chosenMoves.find(nextMove) ==
00345 myitstate->chosenMoves.end()) {
00346 myitstate->unchosenValidMoves.push_back(nextMove);
00347 }
00348 nextMove = origin->nextValidSingleMove(nextMove);
00349 }
00350 }
00351
00352 State*
00353 S_RNAfe_SingleM::RandomNeighborList::first(
00354 State::NeighborList::ItState** itstate) const
00355 {
00356 if (*itstate != NULL) delete *itstate;
00357
00358 *itstate = new ItState();
00359 S_RNAfe_SingleM* rnaState =
00360 new S_RNAfe_SingleM(*origin);
00361
00362 return next(*itstate, rnaState);
00363 }
00364
00365 State*
00366 S_RNAfe_SingleM::RandomNeighborList::next(
00367 State::NeighborList::ItState* itstate, State* elem) const
00368 {
00369 assertbiu(elem != NULL, "Elem is not allowed to be NULL");
00370
00371 ItState* myitstate = static_cast<ItState*>(itstate);
00372 S_RNAfe_SingleM* rnaState =
00373 dynamic_cast<S_RNAfe_SingleM*>(elem);
00374
00375 if (myitstate->switchModeValue > 0) {
00376
00377 size_t pos1, pos2;
00378 SingleMove nextMove;
00379 bool foundValidMove = false;
00380
00381 while (!foundValidMove) {
00382 pos1 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
00383 pos2 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
00384 if (pos1 == pos2) {
00385 continue;
00386 }
00387 nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
00388
00389 if (origin->rnaFreeEnergy.isValidSingleMove(nextMove)) {
00390 if (!(myitstate->chosenMoves.insert(nextMove).second)) {
00391
00392 myitstate->switchModeValue--;
00393 if (myitstate->switchModeValue == 0) {
00394 switchMode(myitstate);
00395 break;
00396 }
00397 }
00398
00399
00400 rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
00401 rnaState->rnaFreeEnergy.applySingleMoveInPlace(
00402 nextMove);
00403 foundValidMove = true;
00404 }
00405 }
00406 }
00407 if (myitstate->switchModeValue == 0) {
00408 size_t n = myitstate->unchosenValidMoves.size();
00409 if (n == 0) {
00410 return NULL;
00411 }
00412 else {
00413 size_t i = biu::RNF::getRN() % n;
00414
00415 rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
00416 rnaState->rnaFreeEnergy.applySingleMoveInPlace(
00417 myitstate->unchosenValidMoves[i]);
00418 myitstate->unchosenValidMoves[i] =
00419 myitstate->unchosenValidMoves[n-1];
00420 myitstate->unchosenValidMoves.resize(n-1);
00421 }
00422 }
00423 return rnaState;
00424 }
00425
00426 State*
00427 S_RNAfe_SingleM::getRandomNeighbor(State* inPlaceNeigh) const
00428 {
00429 S_RNAfe_SingleM* neigh = NULL;
00430
00431 if (inPlaceNeigh != NULL) {
00432 neigh = dynamic_cast<S_RNAfe_SingleM*>(inPlaceNeigh);
00433 assertbiu(neigh != NULL, "inPlaceNeigh is no S_RNAfe_SingleM");
00434 (*neigh) = (*this);
00435 } else {
00436 neigh = new S_RNAfe_SingleM(*this);
00437 }
00438
00439 size_t pos1, pos2;
00440 SingleMove nextMove;
00441
00442 while (true) {
00443 pos1 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
00444 pos2 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
00445 if (pos1 == pos2) {
00446 continue;
00447 }
00448 nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
00449
00450 if (rnaFreeEnergy.isValidSingleMove(nextMove)) {
00451
00452 neigh->rnaFreeEnergy.applySingleMoveInPlace( nextMove );
00453 break;
00454 }
00455 }
00456 return neigh;
00457 }
00458
00459
00460
00462
00463
00464
00465
00466 CSequence
00467 S_RNAfe_SingleM::compress(void) const {
00468 return rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
00469 }
00470
00471
00472
00473
00474
00475
00476 CSequence&
00477 S_RNAfe_SingleM::compress(CSequence& toFill) const {
00478 toFill = rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
00479 return toFill;
00480 }
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490 State*
00491 S_RNAfe_SingleM::uncompress(const CSequence& cseq, State* toFill) const {
00492
00493 biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureRef().size());
00494
00495 if (str.size() != rnaFreeEnergy.getStructureRef().size())
00496 return NULL;
00497
00498 S_RNAfe_SingleM* ret = NULL;
00499
00500 if (toFill == NULL) {
00501 ret = this->clone();
00502 } else {
00503 ret = dynamic_cast<S_RNAfe_SingleM*>(toFill);
00504 assertbiu(ret!=NULL, "given State toFill is no S_RNAfe_SingleM instance");
00505 }
00506
00507 ret->rnaFreeEnergy.setStructure(str);
00508
00509 return ret;
00510 }
00511
00512
00513
00514
00515
00516
00517 State*
00518 S_RNAfe_SingleM::uncompress(const CSequence& cseq) {
00519
00520 biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureRef().size());
00521
00522 if (str.size() != rnaFreeEnergy.getStructureRef().size())
00523 return NULL;
00524
00525 rnaFreeEnergy.setStructure(str);
00526
00527 return this;
00528 }
00529
00530
00531 }