Охереть как много чего сделал

diff --git a/src/dssp.cpp b/src/dssp.cpp
index f28ec01b9cad757560dd0a73e5b9864a76891304..91c562e4d225c25480d55eefc4c0548dfe4506ab 100644 (file)
--- a/src/dssp.cpp
+++ b/src/dssp.cpp
@@ -75,12 +75,9 @@ public:
    void writeOutput() override;
 
 private:
-   real                                                 cutoff_ = 4.0;
-   modeType                                             _modeType;
-   Selection                                            sel_;
+   initParameters                                       initParams;
    std::string                                          filename_;
    int                                                  nres;
-   bool                                                 alternateSearch, verbose;
    DsspTool                                             DT;
    std::vector<std::pair<int, std::string>>             data;
 };
@@ -95,17 +92,17 @@ void Dssp::initOptions(IOptionsContainer* options, TrajectoryAnalysisSettings* s
    };
    options->addOption(
             StringOption("o").store(&filename_).required().defaultValue("SSP.dat").description("Filename for DSSP output"));
-   options->addOption(RealOption("cutoff").store(&cutoff_).required().defaultValue(1.0).description(
+   options->addOption(RealOption("cutoff").store(&initParams.cutoff_).required().defaultValue(4.0).description(
            "cutoff for neighbour search"));
    options->addOption(
-             SelectionOption("sel").store(&sel_).required().description(
+             SelectionOption("sel").store(&initParams.sel_).required().description(
                    "Group for DSSP"));
    options->addOption(
-               EnumOption<modeType>("Mode").store(&_modeType).required().defaultValue(modeType::Classique).enumValue(modeTypeNames).description("Mode for DSSP"));
+               EnumOption<modeType>("HBond Energy Calc Mode").store(&initParams.mode).required().defaultValue(modeType::Classique).enumValue(modeTypeNames).description("Mode for DSSP"));
    options->addOption(
-            BooleanOption("alt").store(&alternateSearch).defaultValue(false).description("Use non-NBsearch method"));
+               EnumOption<NBSearchMethod>("NBSearchMethod Algo").store(&initParams.NBS).required().defaultValue(NBSearchMethod::Classique).enumValue(NBSearchMethodNames).description("Mode for DSSP"));
    options->addOption(
-               BooleanOption("v").store(&verbose).defaultValue(false).description(">:("));
+               BooleanOption("v").store(&initParams.verbose).defaultValue(false).description(">:("));
    settings->setHelpText(desc);
 }
 
@@ -115,7 +112,7 @@ void Dssp::optionsFinished(TrajectoryAnalysisSettings* /* settings */) {}
 
 void Dssp::initAnalysis(const TrajectoryAnalysisSettings &settings, const TopologyInformation& top)
 {
-    DT.setInitValues(cutoff_, sel_, _modeType, verbose);
+    DT.setInitValues(initParams);
     DT.initAnalysis(top);
 }
 

diff --git a/src/dssptools.cpp b/src/dssptools.cpp
index 5e02fe6ef9cf0d21b6e9e7f85610b6d3e7c4c9c6..d59aea9fbc5ce7946941acb66b0c0669ee13e75a 100644 (file)
--- a/src/dssptools.cpp
+++ b/src/dssptools.cpp
@@ -72,15 +72,15 @@ void ResInfo::setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex)
    _ResInfo.at(static_cast<std::size_t>(atomTypeName)) = atomIndex;
 }
 void ResInfo::setInfo(const t_resinfo &info){
-    resinfo = &info;
+    *resinfo = info;
 }
 std::size_t ResInfo::getIndex(backboneAtomTypes atomTypeName) const
 {
    return _ResInfo[static_cast<std::size_t>(atomTypeName)];
 }
 
-t_resinfo ResInfo::getInfo() const{
-    return *resinfo;
+t_resinfo* ResInfo::getInfo() const{
+    return resinfo;
 }
 
 void ResInfo::setHBondedPair(const t_resinfo &PairResinfo){
@@ -93,15 +93,14 @@ std::vector<const t_resinfo*>   ResInfo::getHBondedPairs() const{
 
 secondaryStructures::secondaryStructures(){
 }
-
-void secondaryStructures::initiateSearch(const std::vector<std::vector<bool>> &HBondzMap, const bool PiHelicesPreferenz = true){
+void secondaryStructures::initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez){
     SecondaryStructuresStatusMap.resize(0);
     SecondaryStructuresStringLine.resize(0);
     std::vector<std::size_t> temp; temp.resize(0),
-    PiHelixPreference = PiHelicesPreferenz;
-    HBondsMap = &HBondzMap;
-    SecondaryStructuresStatusMap.resize(HBondsMap.front().size());
-    SecondaryStructuresStringLine.resize(HBondsMap.front().size(), '~');
+    PiHelixPreference = PiHelicesPreferencez;
+    ResInfoMap = &ResInfoMatrix;
+    SecondaryStructuresStatusMap.resize(ResInfoMatrix.size());
+    SecondaryStructuresStringLine.resize(ResInfoMatrix.size(), '~');
 }
 
 void secondaryStructures::secondaryStructuresData::setStatus(const secondaryStructureTypes secondaryStructureTypeName, const bool secondaryStructureStatus = true){
@@ -127,7 +126,16 @@ secondaryStructureTypes secondaryStructures::secondaryStructuresData::getStatus(
     return SSData;
 }
 
-bool secondaryStructures::NoChainBreaksBetween(const std::size_t Resi1, const std::size_t Resi2){
+bool secondaryStructures::hasHBondBetween(std::size_t resi1, std::size_t resi2) const{
+    for(const t_resinfo* &i : (*ResInfoMap)[resi1].getHBondedPairs()){
+        if (i == (*ResInfoMap)[resi2].getInfo()){
+            return true;
+        }
+    }
+    return false;
+}
+
+bool secondaryStructures::NoChainBreaksBetween(std::size_t Resi1, std::size_t Resi2) const{
     bool flag {true};
     if (Resi1 < Resi2){
         for(std::size_t i{Resi1}; flag and i < Resi2; ++i ){
@@ -150,17 +158,17 @@ bool secondaryStructures::NoChainBreaksBetween(const std::size_t Resi1, const st
 
 }
 
-bridgeTypes secondaryStructures::calculateBridge(const std::size_t i, const std::size_t j, const std::vector<std::vector<bool> > &HBondsMap){
+bridgeTypes secondaryStructures::calculateBridge(std::size_t i, std::size_t j) const{
 
-    if( i < 1 || j < 1 || i + 1 >= HBondsMap.size() || j + 1 >= HBondsMap.size() ){
+    if( i < 1 || j < 1 || i + 1 >= ResInfoMap->size() || j + 1 >= ResInfoMap->size() ){
         return bridgeTypes::None;
     }
 
     if(NoChainBreaksBetween(i - 1, i + 1) && NoChainBreaksBetween(j - 1, j + 1)){
-        if((HBondsMap[i + 1][j] && HBondsMap[j][i - 1]) || (HBondsMap[j + 1][i] && HBondsMap[i][j - 1]) ){
+        if((hasHBondBetween(i + 1, j) && hasHBondBetween(j, i - 1)) || (hasHBondBetween(j + 1, i) && hasHBondBetween(i, j - 1)) ){
             return bridgeTypes::Bridge;
         }
-        else if((HBondsMap[i + 1][j - 1] && HBondsMap[j + 1][i - 1]) || (HBondsMap[j][i] && HBondsMap[i][j]) ){
+        else if((hasHBondBetween(i + 1, j - 1) && hasHBondBetween(j + 1, i - 1)) || (hasHBondBetween(j, i) && hasHBondBetween(i, j)) ){
             return bridgeTypes::AntiBridge;
         }
     }
@@ -168,12 +176,12 @@ bridgeTypes secondaryStructures::calculateBridge(const std::size_t i, const std:
     return bridgeTypes::None;
 }
 
-void secondaryStructures::analyzeBridgesAndLaddersPatterns(const std::vector<std::vector<bool>> &HBondsMap){
+void secondaryStructures::analyzeBridgesAndLaddersPatterns(){
 
     for(std::size_t i {1}; i + 4 < SecondaryStructuresStatusMap.size(); ++i){
 
         for(std::size_t j {i + 3}; j + 1 < SecondaryStructuresStatusMap.size(); ++j ){
-            bridgeTypes type = calculateBridge(i, j, HBondsMap);
+            bridgeTypes type = calculateBridge(i, j);
         }
     }
 
@@ -242,10 +250,10 @@ void secondaryStructures::analyzeBridgesAndLaddersPatterns(const std::vector<std
 //    }
 }
 
-void secondaryStructures::analyzeTurnsAndHelicesPatterns(const std::vector<std::vector<bool>> &HBondsMap){
+void secondaryStructures::analyzeTurnsAndHelicesPatterns(){
     for(const turnsTypes &i : { turnsTypes::Turn_4, turnsTypes::Turn_3, turnsTypes::Turn_5 }){
         for(std::size_t j {0}; j + (static_cast<std::size_t>(i) + strideFactor) < SecondaryStructuresStatusMap.size(); ++j){
-            if ( HBondsMap[j][j + (static_cast<std::size_t>(i) + strideFactor)] && NoChainBreaksBetween(j, j + (static_cast<std::size_t>(i) + strideFactor)) ){
+            if ( hasHBondBetween(j, j + (static_cast<std::size_t>(i) + strideFactor)) && NoChainBreaksBetween(j, j + (static_cast<std::size_t>(i) + strideFactor)) ){
                 SecondaryStructuresStatusMap[j + (static_cast<std::size_t>(i) + strideFactor)].setStatus(HelixPositions::End, i);
 
                 for (std::size_t k {1}; k < (static_cast<std::size_t>(i) + strideFactor); ++k){
@@ -319,9 +327,9 @@ void secondaryStructures::analyzeTurnsAndHelicesPatterns(const std::vector<std::
 
 std::string secondaryStructures::patternSearch(const std::vector<std::vector<bool>> &HBondsMap){
 
-    analyzeBridgesAndLaddersPatterns(HBondsMap);
-    analyzeTurnsAndHelicesPatterns(HBondsMap);
-    analyzePPHelicesPatterns(HBondsMap);
+    analyzeBridgesAndLaddersPatterns();
+    analyzeTurnsAndHelicesPatterns();
+    analyzePPHelicesPatterns();
 
     /*Write Data*/
 
@@ -494,11 +502,12 @@ DsspTool::DsspStorage DsspTool::Storage;
 DsspTool::DsspTool(){
 }
 
-void DsspTool::setInitValues(const real _cutoff, const Selection _sel, const modeType _modeType, const bool _verbose){
-    cutoff_ = _cutoff;
-    sel_ = _sel;
-    modeType_ = _modeType;
-    verbose_ = _verbose;
+void DsspTool::setInitValues(const initParameters &initParams){
+    cutoff_ = initParams.cutoff_;
+    sel_ = initParams.sel_;
+    modeType_ = initParams.mode;
+    NBSM_ = initParams.NBS;
+    verbose_ = initParams.verbose;
 }
 
 void DsspTool::calculateBends(const t_trxframe &fr, const t_pbc *pbc)
@@ -568,8 +577,10 @@ bool DsspTool::isHbondExist(const ResInfo& resA,
     * Hbond if E < -0.5
     */
 
-   if(*(resA.getInfo().name) == "PRO"){
-       std::cerr << "PRO detected! This is Resi№ " << *(resA.getInfo().name) << std::endl;
+   std::string proLINE {*(resA.getInfo()->name)};
+
+   if( proLINE.compare("PRO")){
+       std::cerr << "PRO detected! This is Resi№ " << *(resA.getInfo()->name) << std::endl;
        return false;
    }
    const float kCouplingConstant = 27.888;
@@ -726,7 +737,7 @@ void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, const bo
 
        while (as_.findNextPair()){
            if(isHbondExist(IndexMap[as_.getResiI()], IndexMap[as_.getResiJ()], fr, pbc)){
-               IndexMap[as_.getResiI()].setHBondedPair(IndexMap[as_.getResiJ()].getInfo());
+               IndexMap[as_.getResiI()].setHBondedPair(*(IndexMap[as_.getResiJ()].getInfo()));
            }
        }
    }
@@ -748,12 +759,13 @@ void DsspTool::analyzeFrame(int frnr, const t_trxframe &fr, t_pbc *pbc, const bo
        while (pairSearch.findNextPair(&pair))
        {
            if(isHbondExist(IndexMap[pair.refIndex()], IndexMap[pair.testIndex()], fr, pbc)){
-               IndexMap[pair.refIndex()].setHBondedPair(IndexMap[pair.testIndex()].getInfo());
+               IndexMap[pair.refIndex()].setHBondedPair(*(IndexMap[pair.testIndex()].getInfo()));
            }
        }
+
    }
 
-   PatternSearch.initiateSearch(HBondsMap);
+   PatternSearch.initiateSearch(IndexMap, true); // TODO
    calculateBends(fr, pbc);
 
    Storage.storageData(frnr, PatternSearch.calculateLinuxSucks(HBondsMap));
@@ -767,4 +779,3 @@ std::vector<std::pair<int, std::string>> DsspTool::getData(){
 } // namespace analysismodules
 
 } // namespace gmx
-

diff --git a/src/dssptools.h b/src/dssptools.h
index 2507502f734a6a2e1940f484d0826e2fa865790c..f52688e2f8859cb51e5e5beb2de5366c6611cefd 100644 (file)
--- a/src/dssptools.h
+++ b/src/dssptools.h
@@ -72,7 +72,26 @@ const gmx::EnumerationArray<modeType, const char*> modeTypeNames = {
     { "Classique", "Experimental", "Extra" }
 };
 
-//! BackBone atom types. Should be Ca -> C -> O -> N -> H OR pizdec!
+enum class NBSearchMethod : std::size_t
+{
+    Classique = 0,
+    Experimental,
+    Stoopid,
+    Count
+};
+
+const gmx::EnumerationArray<NBSearchMethod, const char*> NBSearchMethodNames = {
+    { "Classique", "Experimental", "Stoopid" }
+};
+
+struct initParameters {
+    Selection                                            sel_;
+    real                                                 cutoff_; // = 4.0; ???
+    modeType                                             mode;
+    NBSearchMethod                                       NBS;
+    bool                                                 verbose;
+};
+
 enum class backboneAtomTypes : std::size_t
 {
    AtomCA,
@@ -94,13 +113,13 @@ public:
    void   setIndex(backboneAtomTypes atomTypeName, std::size_t atomIndex);
    void   setInfo(const t_resinfo &info);
    std::size_t getIndex(backboneAtomTypes atomTypeName) const;
-   t_resinfo   getInfo() const;
+   t_resinfo*   getInfo() const;
    void   setHBondedPair(const t_resinfo &PairResinfo);
    std::vector<const t_resinfo*>   getHBondedPairs() const;
 
 private:
    std::array<std::size_t, static_cast<std::size_t>(backboneAtomTypes::Count)>  _ResInfo{ 0, 0, 0, 0, 0 };
-   const t_resinfo                                                              *resinfo;
+   t_resinfo                                                              *resinfo;
    std::vector<const t_resinfo*>                                                      HBondedResidues;
 };
 
@@ -145,7 +164,7 @@ enum class bridgeTypes : std::size_t {
 class secondaryStructures{
 public:
     secondaryStructures();
-    void                    initiateSearch(const std::vector<std::vector<bool>> &HBondzMap, const bool PiHelicesPreference);
+    void                    initiateSearch(const std::vector<ResInfo> &ResInfoMatrix, const bool PiHelicesPreferencez = true);
     std::string             patternSearch(const std::vector<std::vector<bool>> &HBondsMap);
 //    void                    setStatus(const std::size_t resiIndex, const secondaryStructureTypes secondaryStructureTypeName, const bool secondaryStructureStatus);
 //    bool                    getStatus(const std::size_t resiIndex, const secondaryStructureTypes secondaryStructureTypeName) const;
@@ -168,7 +187,7 @@ public:
     std::vector<secondaryStructuresData>     SecondaryStructuresStatusMap;
 private:
 
-    const std::vector<std::vector<bool>>                                   *HBondsMap;
+    const std::vector<ResInfo>                                       *ResInfoMap;
 
     const gmx::EnumerationArray<secondaryStructureTypes, const char> secondaryStructureTypeNames = {
        { '=', 'S', 'T', 'I', 'G', 'E', 'B', 'H'} // TODO
@@ -177,10 +196,11 @@ private:
     std::string     SecondaryStructuresStringLine;
     const std::size_t     strideFactor{3};
 
-    bridgeTypes     calculateBridge(const std::size_t resi1, const std::size_t resi2, const std::vector<std::vector<bool>> &HBondsMap);
+    bool            hasHBondBetween(std::size_t resi1, std::size_t resi2) const;
 
-    bool            NoChainBreaksBetween(const std::size_t ResiStart, const std::size_t ResiEnd), isLoop(const std::size_t resiIndex) const, PiHelixPreference;
-    void            analyzeBridgesAndLaddersPatterns(const std::vector<std::vector<bool>> &HBondsMap), analyzeTurnsAndHelicesPatterns(const std::vector<std::vector<bool>> &HBondsMap), analyzePPHelicesPatterns(const std::vector<std::vector<bool>> &HBondsMap);
+    bool            NoChainBreaksBetween(std::size_t ResiStart, std::size_t ResiEnd) const, isLoop(const std::size_t resiIndex) const, PiHelixPreference;
+    bridgeTypes     calculateBridge(std::size_t i, std::size_t j) const;
+    void            analyzeBridgesAndLaddersPatterns(), analyzeTurnsAndHelicesPatterns(), analyzePPHelicesPatterns();
 };
 
 class alternateNeighborhoodSearch{
@@ -211,7 +231,7 @@ class DsspTool
 {
 public:
    DsspTool();
-   void setInitValues(const real _cutoff , const Selection _sel, const modeType _modeType, const bool _verbose);
+   void setInitValues(const initParameters &initParams);
    void initAnalysis(/*const TrajectoryAnalysisSettings &settings,*/const TopologyInformation& top);
 
    void analyzeFrame(int frnr, const t_trxframe& fr, t_pbc* pbc, const bool &alternateSearch);
@@ -221,6 +241,7 @@ private:
    real                                  cutoff_ = 4.0;
    Selection                             sel_;
    modeType                              modeType_;
+   NBSearchMethod                        NBSM_;
    bool                                  verbose_;
    AtomsDataPtr                          atoms_;
    std::string                           filename_;