src/cpsp/gecode/GC_HPThreading.cc

Go to the documentation of this file.
00001 /*
00002  *  Main authors:
00003  *     Martin Mann http://www.bioinf.uni-freiburg.de/~mmann/
00004  *
00005  *  Contributing authors:
00006  *     Sebastian Will http://www.bioinf.uni-freiburg.de/~will/
00007  *
00008  *  Copyright:
00009  *     Martin Mann, 2007
00010  *
00011  *  This file is part of the CPSP-tools package:
00012  *     http://www.bioinf.uni-freiburg.de/sw/cpsp/
00013  *
00014  *  See the file "LICENSE" for information on usage and
00015  *  redistribution of this file, and for a
00016  *     DISCLAIMER OF ALL WARRANTIES.
00017  *
00018  */
00019 
00020 
00021 
00022 #include "GC_HPThreading.hh"
00023 #include <biu/Timer.hh>
00024 #include <biu/util/Util_String.h>
00025 #include "cpsp/HCoreDatabase.hh"
00026 #include "GC_ThreadingSpace.hh"
00027 
00028 #include <iostream>
00029 #include <sstream>
00030 #include <iomanip>
00031 
00032 namespace cpsp
00033 {
00034   namespace gecode
00035   {
00036     GC_HPThreading::GC_HPThreading(const biu::LatticeDescriptor* latDescr_,
00037             HCoreDatabase* db, const std::string seq) throw(cpsp::Exception)
00038       :
00039         verboseOutput(false),
00040         latDescr(NULL),
00041         latFrame(NULL),
00042         coreDB(NULL),
00043         sequence(""),
00044         hullLvl(),
00045         coreSize(0),
00046         minEvenOddHs(0),
00047         contactsInSeq(0),
00048         maxPCoreDistance(0),
00049         seqFeatureMap(),
00050         useSymmetryBreaking(true),
00051         useDDS(false),
00052         maxStructures(1),
00053         coreSelection(FIRST),
00054         runMode(ENUMERATE),
00055         noOutput(false),
00056         statisticsOut(false),
00057         noRecomputation(false),
00058         distanceMode(DIST_NONE),
00059         normalizeStructures(false)
00060     {
00061         setLatticeDescriptor(latDescr_);
00062         setSequence(seq);
00063         setCoreDB(db);
00064     }
00065 
00066     GC_HPThreading::~GC_HPThreading()
00067     {
00068         if(latFrame != NULL) { delete latFrame; latFrame = NULL; }
00069 
00070     }
00071 
00072     void
00073     GC_HPThreading::setLatticeDescriptor( const biu::LatticeDescriptor* val )
00074     {
00075         latDescr = val;
00076         if(latFrame != NULL) { delete latFrame; latFrame = NULL; }
00077         if (latDescr != NULL)
00078             latFrame = new biu::LatticeFrame(latDescr, 1);
00079     }
00080 
00081     void
00082     GC_HPThreading::setSequence( const std::string& seq ) throw(cpsp::Exception)
00083     {
00084         sequence = biu::util::Util_String::str2upperCase(seq);
00085         coreSize = 0;
00086         minEvenOddHs = 0;
00087         contactsInSeq = 0;
00088         maxPCoreDistance = 0;
00089 
00090         if (sequence.size() == 0)
00091             return;
00092             // gueltigkeit pruefen
00093         biu::Alphabet hpAlph("HP",1);
00094 //      if (!hpAlph.isAlphabetString(sequence) || sequence.size() < 3) {
00095         if (!hpAlph.isAlphabetString(sequence)) {
00096             throw cpsp::Exception(
00097                 "Error in arguments: sequence is no valid HP-sequence",-2);
00098         }
00099 
00100         size_t i=0, pHead=0, pTail=0, pSubSeq=0, actLvl=0;
00101 
00102             // p-kopf und p-schwanz bestimmen
00103         pHead = sequence.find_first_not_of('P');
00104         if (pHead == std::string::npos) pHead = 0;
00105         pTail = sequence.find_last_not_of('P');
00106         pTail = (pTail == std::string::npos || pTail <= pHead)? 0 : sequence.size() - pTail;
00107 
00108             // core size init
00109         coreSize=(pHead<sequence.size())?1:0;
00110         if (coreSize == 1 && pHead%2==0)    // dann erste H-pos gerade
00111             minEvenOddHs = 1;
00112         else
00113             minEvenOddHs = 0;
00114             // huellen level
00115         hullLvl = GC_ThreadingSpace::HullLevel(sequence.size(),0);  // init
00116             // p-kopf setzen
00117         for (i=0; i<pHead; i++) {
00118             hullLvl[i] = pHead-i;
00119         }
00120             // erster hullLvl init
00121         for (i=pHead+1; i<sequence.size(); i++) {
00122             if (sequence.at(i) == 'H') {                // ende von p-subseq
00123                 pSubSeq = std::max(pSubSeq, actLvl);    // p-subseq setzen
00124                 actLvl = 0;                             // reinit
00125                 coreSize++;                             // kerngroesse ++
00126                 if (i%2==0) minEvenOddHs++;             // gerade H position
00127             } else {
00128                 actLvl++;
00129             }
00130             hullLvl[i] = actLvl;            // level setzen
00131         }
00132             // rueckwaerts angleichen
00133         for (i=sequence.size()-pTail-1; i>pHead+1; i--) {
00134             if (hullLvl[i-1] > (hullLvl[i]+1)) {
00135                 hullLvl[i-1] = hullLvl[i]+1;
00136             }
00137         }
00138             // min von even und odd H positionen bestimmen
00139         minEvenOddHs = std::min( minEvenOddHs, coreSize-minEvenOddHs);
00140 
00141             // singlet positionen bestimmen
00142             // P-singlets
00143         i = sequence.find("PHP");
00144         seqFeatureMap[GC_ThreadingSpace::H_SINGLET].clear();
00145         while ( i != std::string::npos ) {
00146             seqFeatureMap[GC_ThreadingSpace::H_SINGLET].push_back(i+1);
00147             i = sequence.find("PHP",i+2);
00148         }
00149             // h-singlets
00150         i = sequence.find("HPH");
00151         seqFeatureMap[GC_ThreadingSpace::P_SINGLET].clear();
00152         while ( i !=  std::string::npos ) {
00153             seqFeatureMap[GC_ThreadingSpace::P_SINGLET].push_back(i+1);
00154             i = sequence.find("HPH",i+2);
00155         }
00156 
00157             // sequence internal HH-contacts
00158         contactsInSeq = 0;
00159         i = sequence.find("HH");
00160         while ( i !=  std::string::npos ) {
00161             contactsInSeq++;
00162             i = sequence.find("HH", i+1);
00163         }
00164 
00165         maxPCoreDistance = 2* ( std::max(   (pSubSeq/2) +1,
00166                                             std::max( pTail, pHead))
00167                                 +2);
00168 
00169     }
00170     
00171 
00173     unsigned int
00174     GC_HPThreading::enumerateStructures(SuperSpace * space,
00175         unsigned int maxNum, Gecode::Search::Statistics *stat, 
00176         std::ostream &out)
00177     {
00178             // no enumeration using DDS currently
00179         if (useDDS) 
00180             return 0;
00181             // solution number
00182         unsigned int sols = 0;
00183         switch (distanceMode) {
00184             case DIST_NONE : {              
00185                     GC_ThreadingSpace* s = dynamic_cast<GC_ThreadingSpace*>(space);
00186                     assert( s != NULL );
00187                     return enumerateStructures<Gecode::DFS<GC_ThreadingSpace>,GC_ThreadingSpace> 
00188                         ( s, maxNum, stat, out); 
00189                     break;
00190                 }
00191             case DIST_SHAPES : {
00192                     GC_ThreadingSpaceShapes* s = dynamic_cast<GC_ThreadingSpaceShapes*>(space);
00193                     assert( s != NULL );
00194                     return enumerateStructures<Gecode::DFS<GC_ThreadingSpaceShapes>,GC_ThreadingSpaceShapes> 
00195                         ( s, maxNum, stat, out); 
00196                     break;
00197                 }
00198             case DIST_POSITIONS : {
00199                     GC_ThreadingSpacePosDist* s = dynamic_cast<GC_ThreadingSpacePosDist*>(space);
00200                     assert( s != NULL );
00201                     return enumerateStructures<Gecode::BAB<GC_ThreadingSpacePosDist>,GC_ThreadingSpacePosDist> 
00202                         ( s, maxNum, stat, out); 
00203                     break;
00204                 }
00205             case DIST_ABSMOVES : {
00206                     GC_ThreadingSpaceMoveDist* s = dynamic_cast<GC_ThreadingSpaceMoveDist*>(space);
00207                     assert( s != NULL );
00208                     return enumerateStructures<Gecode::BAB<GC_ThreadingSpaceMoveDist>,GC_ThreadingSpaceMoveDist> 
00209                         ( s, maxNum, stat, out); 
00210                     break;
00211                 }
00212             default : 
00213                     assertbiu(false,"distance mode not supported");
00214                     break;
00215         }
00216         return sols;
00217     }
00218 
00219     unsigned int
00220     GC_HPThreading::countStructures(SuperSpace * space,
00221         unsigned int maxNum, bool useDDS, Gecode::Search::Statistics *actStat)
00222     {
00223         switch(distanceMode) {
00224 
00225                 // normal counting
00226             case DIST_NONE : {
00227                 GC_ThreadingSpace* s = dynamic_cast<GC_ThreadingSpace*>(space);
00228                 assert(s != NULL);
00229                 if (useDDS) {   // using DDS
00230                     assertbiu(false,"DDS not supported");
00231                 } else {        // using DFS
00232                     return Gecode::countDFS<unsigned int, GC_ThreadingSpace> (
00233                                 s,
00234                                 (noRecomputation?1:Gecode::Search::Config::c_d),
00235                                 (noRecomputation?1:Gecode::Search::Config::a_d),
00236                                 actStat, NULL,
00237                                 maxNum
00238                             );
00239                 }
00240                 break;
00241             }
00242 
00243                 // normal counting
00244             case DIST_SHAPES : {
00245                 GC_ThreadingSpaceShapes* s = dynamic_cast<GC_ThreadingSpaceShapes*>(space);
00246                 assert(s != NULL);
00247                 if (useDDS) {   // using DDS
00248                     assertbiu(false,"DDS not supported");
00249                 } else {        // using DFS
00250                     return Gecode::countDFS<unsigned int, GC_ThreadingSpaceShapes> (
00251                                 s,
00252                                 (noRecomputation?1:Gecode::Search::Config::c_d),
00253                                 (noRecomputation?1:Gecode::Search::Config::a_d),
00254                                 actStat, NULL,
00255                                 maxNum
00256                             );
00257                 }
00258                 break;
00259             }
00260 
00261                 // BAB counting introducing distance constraints to already 
00262                 // found solutions on positions     
00263             case DIST_POSITIONS : {
00264                 GC_ThreadingSpacePosDist* s = dynamic_cast<GC_ThreadingSpacePosDist*>(space);
00265                 assert(s != NULL);
00266                 if (useDDS) {
00267                     assertbiu(false,"DDS not supported");
00268                 } else {
00269                         // using branch and bound
00270                     return countBAB<unsigned int, GC_ThreadingSpacePosDist> ( 
00271                                 s,
00272                                 (noRecomputation?1:Gecode::Search::Config::c_d),
00273                                 (noRecomputation?1:Gecode::Search::Config::a_d),
00274                                 actStat, NULL,
00275                                 maxNum
00276                             );
00277                 } 
00278                 break;
00279             }
00280                 // BAB counting introducing distance constraints to already 
00281                 // found solutions on absolute moves                
00282             case DIST_ABSMOVES : {
00283                 GC_ThreadingSpaceMoveDist* s = dynamic_cast<GC_ThreadingSpaceMoveDist*>(space);
00284                 assert(s != NULL);
00285                 if (useDDS) {
00286                     assertbiu(false,"DDS not supported");
00287                 } else {
00288                         // using branch and bound
00289                     return countBAB<unsigned int, GC_ThreadingSpaceMoveDist> ( 
00290                                 s,
00291                                 (noRecomputation?1:Gecode::Search::Config::c_d),
00292                                 (noRecomputation?1:Gecode::Search::Config::a_d),
00293                                 actStat, NULL,
00294                                 maxNum
00295                             );
00296                 } 
00297                 break;
00298             }
00299             default : { 
00300                     assertbiu(false,"distance mode not supported");
00301                 break;
00302             }
00303         }
00304         return 0;
00305     }
00306 
00307     unsigned int
00308     GC_HPThreading::generateStructures(void) {
00309 
00310         assert(latFrame != NULL);
00311         assert(coreDB != NULL);
00312         
00313         if (sequence.size() == 0)
00314             return 0;
00315         
00316         HCore core;
00317         unsigned int sols = 0, oldSols = 0, oldContacts = 0,
00318             oldContactSols = 0, maxContactsFound = 0;
00319 
00320         if (!noOutput)
00321             std::cout <<"==> HP-sequence = " <<sequence <<std::endl;
00322 
00323         Gecode::Search::Statistics *globStat= NULL, *actStat=NULL;
00324         if (verboseOutput) {
00325             std::cout << "\ngenerating structures via Gecode\n" << std::endl;
00326         }
00327         if (statisticsOut) {
00328             globStat = new Gecode::Search::Statistics;
00329             actStat = new Gecode::Search::Statistics;
00330         }
00331             // solution storage for enumeration with minimal distance
00332         GC_ThreadingSpacePosDist::SolutionVec* solPos = NULL, *solMove = NULL;
00333 
00334             // timer setzen
00335         biu::Timer timer; timer.start();
00336         biu::Timer sTimer;
00337         double searchTime= 0.0;
00338 
00339         GC_ThreadingSymmBreaker::GlobalShiftVec* symmetryShiftVec = NULL;
00340             // neighborhood vectors for neighboring constraints
00341         biu::IndexVec indexedNeighVecs;
00342          
00343             // set maximal frame size needed
00344         latFrame->setFrameSize( maxPCoreDistance
00345                                 + coreDB->getActCoreSize());
00346 
00347         while (oldSols < maxStructures && coreDB->getNextCore(core)) {
00348             try {
00349                     // guide output and energy level control
00350                 if (core.getContacts() > oldContacts) {
00351                     oldContacts = core.getContacts();
00352                     if (!noOutput)
00353                         std::cout   <<"\n==> energy = -"
00354                                     <<(oldContacts - contactsInSeq)
00355                                     <<std::endl;
00356                 } 
00357                 
00358                 if (core.getContacts() < oldContacts) {
00359                     if (!noOutput)
00360                         std::cout   <<"\n ==> number generated structures = "
00361                                     <<(oldSols - oldContactSols)
00362                                     <<std::endl;
00363                         // number of solutions since last energy level change
00364                     oldContactSols = oldSols;
00365                     oldContacts = core.getContacts();
00366                     if (!noOutput)
00367                         std::cout   <<"\n==> energy = -"
00368                                     <<(oldContacts - contactsInSeq)
00369                                     <<std::endl;
00370                 }
00371                     // set minimal frame size
00372                 if (distanceMode == DIST_NONE ) {   
00373                         // not possible if distance of solution is forced 
00374                         // (different positions indices etc.)  
00375                     latFrame->setFrameSize( maxPCoreDistance
00376                                             + core.getMaxDimension());
00377                 }
00378 
00379                 if (verboseOutput) {
00380                     std::cout <<" => next core ("<<core.getContacts()
00381                         <<" contacts, "
00382                         <<core.getNumOfPSinglets(latFrame)
00383                         <<" P-singlets) "
00384                         <<std::endl;
00385                 }
00386                     // naechster, wenn nicht genug singlet-positionen vorhanden
00387                 if ( core.getNumOfPSinglets(latFrame)
00388                             < (int)  seqFeatureMap[GC_ThreadingSpace::P_SINGLET].size()
00389                     || core.getNumOfHSinglets(latFrame)
00390                             < (int)  seqFeatureMap[GC_ThreadingSpace::H_SINGLET].size())
00391                 {
00392                     if (verboseOutput) {
00393                         std::cout   <<"     --> skipped (singlet positions)"
00394                                     <<std::endl;
00395                     }
00396                     continue;
00397                 }
00398                     // naechster, wenn cubisches gitter und falsches verhaeltnis
00399                     // zwischen gerade und ungeraden kernpositionen
00400                 if ( latFrame->getDescriptor()->getName() == "cub"
00401                     && core.getMinEvenOdd() != minEvenOddHs )
00402                 {
00403                     if (verboseOutput) {
00404                         std::cout   <<"     --> skipped (cubic even/odd ratio) "
00405                                     <<std::endl;
00406                     }
00407                     continue;
00408                 }
00409 
00410                     // generate global shift vectors for symmetry breaking
00411                 if (useSymmetryBreaking) {
00412                     symmetryShiftVec = GC_ThreadingSymmBreaker::generateGlobalShiftVec(
00413                                         &core, latFrame, symmetryShiftVec);
00414                 }
00415                 
00416                 indexedNeighVecs.clear();
00417                 indexedNeighVecs = latFrame->getIndexedNeighborhood();
00418 
00419                     // new space
00420                 GC_ThreadingSpace* s = NULL;
00421                 int branching = BR_NONE;
00422                 switch (distanceMode) {
00423                     case DIST_NONE : 
00424                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00425                         s = new GC_ThreadingSpace(&sequence,
00426                         latFrame, &indexedNeighVecs, &seqFeatureMap, &hullLvl, &core,
00427                         symmetryShiftVec, branching);
00428                         break;
00429                     case DIST_SHAPES : 
00430                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00431                         s = new GC_ThreadingSpaceShapes(&sequence,
00432                         latFrame, &indexedNeighVecs, &seqFeatureMap, &hullLvl, &core,
00433                         symmetryShiftVec, branching);
00434                         break;
00435                     case DIST_POSITIONS :
00436                         if (solPos == NULL)
00437                             solPos = new GC_ThreadingSpacePosDist::SolutionVec();
00438                         else
00439                             solPos->clear(); // removed found sequences because core changed 
00440                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00441                         s = new GC_ThreadingSpacePosDist(&sequence,
00442                             latFrame, &indexedNeighVecs, &seqFeatureMap, 
00443                             &hullLvl, &core, symmetryShiftVec, 
00444                             branching,
00445                             solPos, 
00446                             (std::max(sequence.size()-distanceMin,(std::string::size_type)0)) );
00447                         break;
00448                     case DIST_ABSMOVES :
00449                         if (solMove == NULL)
00450                             solMove = new GC_ThreadingSpaceMoveDist::SolutionVec(); 
00451                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00452                         s = new GC_ThreadingSpaceMoveDist(&sequence,
00453                             latFrame, &indexedNeighVecs, &seqFeatureMap, 
00454                             &hullLvl, &core, symmetryShiftVec, 
00455                             branching,
00456                             solMove, 
00457                             (std::max(sequence.size()-1-distanceMin,(std::string::size_type)0)) );
00458                         break;
00459                     default :
00460                         break;
00461                 }
00462 
00463                 sTimer.start();
00464                 if (runMode == ENUMERATE) {
00465                     sols = enumerateStructures(s, maxStructures-oldSols, actStat);
00466                 } else {
00467                     sols = countStructures(s, maxStructures-oldSols, useDDS, actStat);
00468                 }
00469                 searchTime += sTimer.stop();
00470                 if (globStat != NULL) {
00471                     globStat->clone     += actStat->clone;
00472                     globStat->fail      += actStat->fail;
00473                     globStat->commit    += actStat->commit;
00474 //                  globStat->decomp    += actStat->decomp;
00475                     globStat->propagate += actStat->propagate;
00476                     globStat->memory    += actStat->memory;
00477                 }
00478                     // increase maximal found contacts
00479                 if (sols > 0) {
00480                     maxContactsFound = std::max( maxContactsFound, 
00481                                                 (oldContacts - contactsInSeq));
00482                 }
00483 
00484                 oldSols += sols;
00485                 delete s;
00486                     // wenn nur beste gesucht werden, dann u.u. energieschranke
00487                     // anpassen
00488                 if (    coreSelection == ALL_BEST
00489                         && sols > 0) {
00490                     if (!noOutput)
00491                         std::cout   <<"   + " <<(sols)
00492                                     <<" structures" <<std::endl;
00493                     coreDB->setActMinHHcontacts(core.getContacts());
00494                 }
00495                 if (symmetryShiftVec != NULL) { delete symmetryShiftVec; symmetryShiftVec = NULL; }
00496             } catch (Gecode::Exception e) {
00497                 std::cout << "Exception in Gecode: "
00498                  << e.what() << "." << std::endl
00499                 << "Stopping..." << std::endl;
00500             }
00501         }
00502         indexedNeighVecs.clear();
00503         if (symmetryShiftVec != NULL) { delete symmetryShiftVec; symmetryShiftVec = NULL; }
00504         if (!noOutput)
00505             std::cout   <<"\n ==> number generated structures = "
00506                         <<(oldSols - oldContactSols)
00507                         <<"\n\n- total number of generated structures = "
00508                         <<std::setw(15) <<oldSols
00509                         <<"\n- maximal contacts reachable           = "
00510                         <<std::setw(15) << maxContactsFound
00511                         <<std::endl;
00512 
00513         if (globStat != NULL) {
00514                 // statistik ausgabe
00515             std::cout << "- statistics"
00516             << "\n\tsolutions:    " << oldSols
00517             << "\n\tpropagations: " << globStat->propagate 
00518             << "\n\tfailures:     " << globStat->fail
00519             << "\n\tclones:       " << globStat->clone 
00520 //          << "\n\tdecompos.:    " << globStat->decomp 
00521             << "\n\tcommits:      " << globStat->commit 
00522             << "\n\tpeak memory:  "
00523                 << static_cast<int>((globStat->memory+1023) / 1024) << " KB"
00524             << "\n\tsearchtime:   " << searchTime <<" ms" 
00525             << "\n\truntime:      " << timer.stop() <<" ms"
00526             << std::endl;
00527             
00528                 // garbage collection
00529             delete globStat; globStat = NULL;
00530             delete actStat;  actStat = NULL;
00531         }
00532         
00533         if (solPos != NULL) { solPos->clear(); delete solPos; solPos = NULL; }
00534         if (solMove != NULL) { solMove->clear(); delete solMove; solMove = NULL; }
00535 
00536         return oldSols;
00537     }
00538 
00539         // generates structures and stores them in structs
00540     unsigned int
00541     GC_HPThreading::enumerateStructures(std::vector<std::string> & structs) {
00542 
00543         assert(latFrame != NULL);
00544         assert(coreDB != NULL);
00545         
00546         structs.clear();
00547         
00548         if (sequence.size() == 0)
00549             return 0;
00550 
00551         HCore core;
00552         unsigned int sols = 0, oldSols = 0, oldContacts = INT_MAX,
00553             oldContactSols = 0;
00554 
00555         GC_ThreadingSymmBreaker::GlobalShiftVec* symmetryShiftVec = NULL;
00556             // neighborhood vectors for neighboring constraints
00557         biu::IndexVec indexedNeighVecs;
00558 
00559             // solution storage for enumeration with minimal distance
00560         GC_ThreadingSpacePosDist::SolutionVec* solPos = NULL, *solMove = NULL;
00561 
00562             // set maximal frame size needed
00563         latFrame->setFrameSize( maxPCoreDistance
00564                                 + coreDB->getActCoreSize());
00565 
00566         while (oldSols < maxStructures && coreDB->getNextCore(core)) {
00567             try {
00568                     // set minimal frame size
00569                     // set minimal frame size
00570                 if (distanceMode == DIST_NONE ) {   
00571                         // not possible if distance of solution is forced 
00572                         // (different positions indices etc.)  
00573                     latFrame->setFrameSize( maxPCoreDistance
00574                                             + core.getMaxDimension());
00575                 }
00576 
00577                     // naechster, wenn nicht genug singlet-positionen vorhanden
00578                 if ( core.getNumOfPSinglets(latFrame)
00579                             < (int)  seqFeatureMap[GC_ThreadingSpace::P_SINGLET].size()
00580                     || core.getNumOfHSinglets(latFrame)
00581                             < (int)  seqFeatureMap[GC_ThreadingSpace::H_SINGLET].size())
00582                 {
00583                     continue;
00584                 }
00585                     // naechster, wenn cubisches gitter und falsches verhaeltnis
00586                     // zwischen gerade und ungeraden kernpositionen
00587                 if ( latFrame->getDescriptor()->getName() == "cub"
00588                     && core.getMinEvenOdd() != minEvenOddHs )
00589                 {
00590                     continue;
00591                 }
00592 
00593                     // generate global shift vectors for symmetry breaking
00594                 if (useSymmetryBreaking) {
00595                     symmetryShiftVec = GC_ThreadingSymmBreaker::generateGlobalShiftVec(
00596                                         &core, latFrame, symmetryShiftVec);
00597                 }
00598                 
00599                 indexedNeighVecs.clear();
00600                 indexedNeighVecs = latFrame->getIndexedNeighborhood();
00601 
00602                     // new space
00603                     // new space
00604                 SuperSpace* s = NULL;
00605                 int branching = BR_NONE;
00606                 switch (distanceMode) {
00607                     case DIST_NONE : 
00608                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00609                         s = new GC_ThreadingSpace(&sequence,
00610                         latFrame, &indexedNeighVecs, &seqFeatureMap, &hullLvl, &core,
00611                         symmetryShiftVec, branching);
00612                         break;
00613                     case DIST_SHAPES : 
00614                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00615                         s = new GC_ThreadingSpaceShapes(&sequence,
00616                         latFrame, &indexedNeighVecs, &seqFeatureMap, &hullLvl, &core,
00617                         symmetryShiftVec, branching);
00618                         break;
00619                     case DIST_POSITIONS :
00620                         if (solPos == NULL)
00621                             solPos = new GC_ThreadingSpacePosDist::SolutionVec();
00622                         else
00623                             solPos->clear(); // clear found sols due to new core  
00624                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00625                         s = new GC_ThreadingSpacePosDist(&sequence,
00626                             latFrame, &indexedNeighVecs, &seqFeatureMap, 
00627                             &hullLvl, &core, symmetryShiftVec, 
00628                             branching,
00629                             solPos, 
00630                             (std::max(sequence.size()-distanceMin,(std::string::size_type)0)) );
00631                         break;
00632                     case DIST_ABSMOVES :
00633                         if (solMove == NULL)
00634                             solMove = new GC_ThreadingSpacePosDist::SolutionVec(); 
00635                         branching = (useDDS?BR_DDS:BR_DFS)|(useSymmetryBreaking?BR_SYM:BR_NONE);
00636                         s = new GC_ThreadingSpaceMoveDist(&sequence,
00637                             latFrame, &indexedNeighVecs, &seqFeatureMap, 
00638                             &hullLvl, &core, symmetryShiftVec, 
00639                             branching,
00640                             solMove, 
00641                             (std::max(sequence.size()-1-distanceMin,(std::string::size_type)0)) );
00642                         break;
00643                     default :
00644                         break;
00645                 }
00646                     // ausgabe der energieniveaus
00647                 if (core.getContacts() < oldContacts) {
00648                         // number of solutions since last energy level change
00649                     oldContactSols = oldSols;
00650                     oldContacts = core.getContacts();
00651                 }
00652                 std::ostringstream str;
00653                 sols = enumerateStructures(s, maxStructures-oldSols, NULL, str);
00654                 std::string retStr = str.str();
00655                 size_t b=0, e=retStr.find('\n');
00656                 do {
00657                     structs.push_back(retStr.substr(b,(e==retStr.npos?retStr.size()-b:e-b)));
00658                     b = e+1;
00659                     e = retStr.find('\n',b);
00660                 } while(e != retStr.npos); 
00661                 
00662                 oldSols += sols;
00663                 delete s;
00664                     // wenn nur beste gesucht werden, dann u.u. energieschranke
00665                     // anpassen
00666                 if (    coreSelection == ALL_BEST && sols > 0) {
00667                     coreDB->setActMinHHcontacts(core.getContacts());
00668                 }
00669                 if (symmetryShiftVec != NULL) { delete symmetryShiftVec; symmetryShiftVec = NULL; }
00670             } catch (Gecode::Exception e) {
00671                 std::cout << "Exception in Gecode: "
00672                  << e.what() << "." << std::endl
00673                 << "Stopping..." << std::endl;
00674             }
00675         }
00676         indexedNeighVecs.clear();
00677         if (symmetryShiftVec != NULL) { delete symmetryShiftVec; symmetryShiftVec = NULL; }
00678         if (solPos != NULL) { solPos->clear(); delete solPos; solPos = NULL; }
00679         if (solMove != NULL) { solMove->clear(); delete solMove; solMove = NULL; }
00680 
00681         return oldSols;
00682     }
00683 
00684 
00685   } // namespace gecode
00686 } // namespace cpsp
00687