src/cpsp/HPThreadingOptions.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 #include "cpsp/HPThreadingOptions.hh"
00021 
00022 #include <limits.h>
00023 
00024 #include <biu/util/Util_String.h>
00025 
00026 #include <biu/LatticeDescriptorCUB.hh>
00027 #include <biu/LatticeDescriptorFCC.hh>
00028 
00029 #include "cpsp/HCoreDatabaseFILE.hh"
00030 //#include "cpsp/HCoreDatabasePGSQL.hh"
00031 
00032 
00033 namespace cpsp
00034 {
00035 
00036     HPThreadingOptions::HPThreadingOptions()  :
00037         sequence(NULL), latDescr(NULL), lattice(NULL), 
00038         coreSelection(FIRST), minE(INT_MIN/2+2), maxE(0), maxSols(0),
00039         hullLvl(), maxPCoreDistance(0), coreDB(NULL),
00040         onlyCountingStructures(false), verboseOutput(false),
00041         minEvenOddHs(0),
00042         contactsInSeq(0),
00043         useSymmetryBreaking(true),
00044         useNoDDS(false),
00045         normalizeStructures(false)
00046         
00047     {
00048     }
00049 
00050     HPThreadingOptions::HPThreadingOptions(int argc, char** argv) 
00051         throw(cpsp::Exception) :
00052         sequence(NULL), latDescr(NULL), lattice(NULL), 
00053         coreSelection(FIRST), minE(INT_MIN/2+2), maxE(0), maxSols(0),
00054         hullLvl(), maxPCoreDistance(0), coreDB(NULL),
00055         onlyCountingStructures(false), verboseOutput(false),
00056         minEvenOddHs(0),
00057         contactsInSeq(0),
00058         useSymmetryBreaking(true),
00059         useNoDDS(false),
00060         normalizeStructures(false)
00061         
00062     {
00063             // init allowed arguments
00064         biu::OptionMap allowedArgs;
00065         std::string infoText;
00066         initAllowedArguments(allowedArgs,infoText);
00067         
00068             // parse programm arguments 
00069         biu::COptionParser opts = biu::COptionParser(   allowedArgs, argc, 
00070                                                         (char**)argv, infoText);
00071         
00072         if (opts.noErrors()) {
00073     
00074             if (opts.getBoolVal("help")) {
00075                 throw cpsp::Exception(
00076                     "",0);
00077             }
00078             
00079                 // mode flags 
00080             verboseOutput = opts.getBoolVal("v");
00081             onlyCountingStructures = opts.getBoolVal("count");
00082             
00083                 // sequence
00084             std::string tmpSeq = opts.getStrVal("seq");
00085             unsigned int coreSize = setSequence(tmpSeq);
00086                         
00087                 // lattice modell
00088             if (opts.getStrVal("lat") == "CUB")
00089                 latDescr = new biu::LatticeDescriptorCUB();
00090             else if (opts.getStrVal("lat") == "FCC")
00091                 latDescr = new biu::LatticeDescriptorFCC();
00092             else 
00093                 throw cpsp::Exception(
00094                     "Error in arguments: lattice is not supported",-2);
00095             
00096                 
00097             lattice = new biu::LatticeFrame(latDescr, 2*sequence->size()+4);
00098             
00099                 // core selection mode
00100                 // <all> overwrites <allbest> if both given
00101             if (opts.getBoolVal("all")) {
00102                 coreSelection = ALL;
00103             } else if (opts.getBoolVal("allbest")) {
00104                 coreSelection = ALL_BEST;
00105             }
00106             if (coreSelection == FIRST) {
00107                 maxSols = 1;
00108             }
00109                 // energie bounds
00110 // @TODO option wieder reinnehmen               
00111 //          minE = opts.getIntVal("minE");
00112 //          maxE = opts.getIntVal("maxE");
00113             
00114                 // contact bounds for database access
00115             unsigned int minContacts = -maxE + contactsInSeq;
00116             unsigned int maxContacts = -minE + contactsInSeq;
00117             
00118                 // check for semantic errors
00119             if (minE > maxE)
00120                 throw cpsp::Exception("Error in arguments : minE > maxE", -2);
00121                 
00122                 // maximal structure number to enumerate
00123             if (coreSelection != FIRST) {
00124                 int tmp = opts.getIntVal("maxSol");
00125                 if (tmp < 1)
00126                     throw cpsp::Exception("Error in arguments : maxSol <= 0", -2);
00127                 maxSols = (unsigned int) tmp;
00128             }
00129 
00130             // use symmetry breaking
00131             useSymmetryBreaking = ! opts.getBoolVal("noSymmBreak");
00132             
00133             // use no dynamic decomposition search (DDS)
00134             // use DDS only in counting mode
00135             useNoDDS = !onlyCountingStructures || opts.getBoolVal("noDDS");
00136             
00137             
00138             normalizeStructures = opts.getBoolVal("normalize");
00139 
00140                 // H-core database
00141 
00142             if (opts.argExist("dbPath")) {
00143                 coreDB = new HCoreDatabaseFILE(opts.getStrVal("dbPath"));
00144             } else {
00145                 throw cpsp::Exception(
00146                     "H-core database Error: no database available (FILE or SQL)",-1);
00147             }
00148                 // test connection
00149             if (!coreDB->isConnected()) {
00150                 throw cpsp::Exception(
00151                     "H-core database Error: can't open database",-1);
00152             }
00153 /*
00154             std::cout <<"coreSize = " <<coreSize <<std::endl;
00155             std::cout <<"minContacts = " <<minContacts <<std::endl;
00156             std::cout <<"maxContacts = " <<maxContacts <<std::endl;
00157             std::cout <<"latDescr = " <<latDescr->getName() <<std::endl;
00158 */          
00159             if (!coreDB->initCoreAccess(    *latDescr,
00160                                             coreSize,
00161                                             minContacts,
00162                                             maxContacts)) {
00163                 throw cpsp::Exception(
00164                     "H-core database Error: access init failed",-1);
00165             }
00166         } else {
00167             throw cpsp::Exception(
00168                 "Parameter Error: not parseable",-1);
00169         }
00170     }
00171     
00172     HPThreadingOptions::~HPThreadingOptions() 
00173     {
00174             // aufraeumen
00175         if( coreDB != NULL) {
00176             delete coreDB; coreDB = NULL;
00177         }
00178         if( lattice != NULL ) {
00179             delete lattice; lattice = NULL;
00180         }
00181         if( latDescr != NULL ) {
00182             delete latDescr; latDescr = NULL;
00183         }
00184         if( sequence != NULL ) {
00185             delete sequence; sequence = NULL;
00186         }
00187     }
00188     
00189     unsigned int 
00190     HPThreadingOptions::setSequence(const std::string &seq) throw(cpsp::Exception)
00191     {
00192         sequence = new std::string(biu::util::Util_String::str2upperCase(seq));
00193             // check if valid
00194         biu::Alphabet hpAlph("HP",1);
00195         if (!hpAlph.isAlphabetString(*sequence) || sequence->size() < 3) {
00196             throw cpsp::Exception(
00197                 "Error in arguments: sequence is no valid HP-sequence",-2);
00198         }
00199         
00200         size_t i=0, pHead=0, pTail=0, pSubSeq=0, actLvl=0;
00201 
00202             // get P-head und P-tail
00203         pHead = sequence->find_first_not_of('P');
00204         pTail = sequence->find_last_of('P');
00205         pTail = (pTail == std::string::npos)? 0 : sequence->size() - pTail;
00206         
00207             // core size init
00208         unsigned int coreSize=(pHead<sequence->size())?1:0;
00209         if (coreSize == 1 && pHead%2==0)    // first H-pos is even
00210             minEvenOddHs = 1;
00211             // hull level
00212         hullLvl = std::vector<unsigned int>(sequence->size(),0);    // init
00213             // set P-head
00214         for (i=0; i<pHead; i++) {
00215             hullLvl[i] = pHead-i;
00216         }
00217             // init first hullLvl 
00218         for (i=pHead+1; i<sequence->size(); i++) {
00219             if (sequence->at(i) == 'H') {               // end of p-subseq
00220                 pSubSeq = std::max(pSubSeq, actLvl);    // set p-subseq 
00221                 actLvl = 0;                             // reinit
00222                 coreSize++;                             // kern size ++
00223                 if (i%2==0) minEvenOddHs++;             // even H position
00224             } else {
00225                 actLvl++;
00226             }
00227             hullLvl[i] = actLvl;            // set hull level 
00228         }
00229             // adjust backwards
00230         for (i=sequence->size()-pTail-1; i>pHead+1; i--) {
00231             if (hullLvl[i-1] > (hullLvl[i]+1)) {
00232                 hullLvl[i-1] = hullLvl[i]+1;
00233             }
00234         }
00235             // get min of even and odd H positions
00236         minEvenOddHs = std::min( minEvenOddHs, coreSize-minEvenOddHs);
00237         
00238             // get singlet position
00239             // P-singlets
00240         i = sequence->find("PHP");
00241         seqFeatureMap[H_SINGLET].clear();
00242         while ( i != std::string::npos ) {
00243             seqFeatureMap[H_SINGLET].push_back(i+1);
00244             i = sequence->find("PHP",i+2);
00245         }
00246             // h-singlets
00247         i = sequence->find("HPH");
00248         seqFeatureMap[P_SINGLET].clear();
00249         while ( i !=  std::string::npos ) {
00250             seqFeatureMap[P_SINGLET].push_back(i+1);
00251             i = sequence->find("HPH",i+2);
00252         }
00253         
00254             // sequence internal HH-contacts
00255         contactsInSeq = 0;
00256         i = sequence->find("HH");
00257         while ( i !=  std::string::npos ) {
00258             contactsInSeq++;
00259             i = sequence->find("HH", i+1);
00260         }
00261         
00262             // set lattice frame size 
00263         maxPCoreDistance = 2* ( std::max(   (pSubSeq/2) +1, 
00264                                             std::max( pTail, pHead))
00265                                 +2);
00266         
00267         return coreSize;
00268     }
00269     
00270     
00271     void
00272     HPThreadingOptions::initAllowedArguments(biu::OptionMap & allowedArgs,
00273                                             std::string &infoText ) const {
00274         allowedArgs.push_back(biu::COption( 
00275                                 "dbPath", false, biu::COption::STRING, 
00276                                 "file based database root path"));
00277         allowedArgs.push_back(biu::COption( 
00278                                 "seq", false, biu::COption::STRING, 
00279                                 "HP-sequence to fold"));
00280         allowedArgs.push_back(biu::COption( 
00281                                 "lat", true, biu::COption::STRING, 
00282                                 "lattice (CUB, FCC)", "CUB"));
00283         allowedArgs.push_back(biu::COption( 
00284                                 "all", true, biu::COption::BOOL, 
00285                                 "search all solutions"));
00286         allowedArgs.push_back(biu::COption( 
00287                                 "noSymmBreak", true, biu::COption::BOOL, 
00288                                 "do not break symmetries"));
00289         allowedArgs.push_back(biu::COption( 
00290                                 "noDDS", true, biu::COption::BOOL, 
00291                                 "do not use Dynamic Decomposition Search"));
00292         allowedArgs.push_back(biu::COption( 
00293                                 "allbest", true, biu::COption::BOOL, 
00294                                 "search all solutions with best energy"));
00295 /*      allowedArgs.push_back(biu::COption( 
00296                                 "minE", true, biu::COption::INT, 
00297                     "lower energy bound for the search (value <= maxE <= 0)",
00298                                 biu::util::Util_String::int2str(INT_MIN/2+2)));
00299         allowedArgs.push_back(biu::COption( 
00300                                 "maxE", true, biu::COption::INT, 
00301                             "upper energy bound for the search (value <= 0)",
00302                                 "0"));
00303 */      allowedArgs.push_back(biu::COption( 
00304                                 "maxSol", true, biu::COption::INT, 
00305                             "maximum of solutions to search for (value > 0)",
00306                                 biu::util::Util_String::int2str(INT_MAX-2)));
00307         allowedArgs.push_back(biu::COption( 
00308                                 "count", true, biu::COption::BOOL, 
00309                                 "only counting structures"));
00310         allowedArgs.push_back(biu::COption( 
00311                                 "normalize", true, biu::COption::BOOL, 
00312                                 "print structures as normalized move strings"));
00313         allowedArgs.push_back(biu::COption( 
00314                                 "v", true, biu::COption::BOOL, 
00315                                 "verbose output"));
00316         allowedArgs.push_back(biu::COption( 
00317                                 "help", true, biu::COption::BOOL, 
00318                                 "program parameters and help"));
00319 
00320         infoText = " - the option <all> overwrites the <allbest> option";
00321     
00322     } // initArguments
00323 
00324 } // namespace cpsp