From 6bf3a6108ea96c468ed9a35728d487f6460a4e41 Mon Sep 17 00:00:00 2001
From: Valerio Mariani <valerio.mariani@unibas.ch>
Date: Wed, 21 Mar 2012 17:38:23 +0100
Subject: [PATCH] Updated LDT to use now distance list implementation

---
 modules/mol/alg/pymod/wrap_mol_alg.cc  |   4 +-
 modules/mol/alg/src/ldt.cc             |   3 +-
 modules/mol/alg/src/local_dist_test.cc | 325 ++++++++++++++++---------
 modules/mol/alg/src/local_dist_test.hh |  57 ++++-
 4 files changed, 261 insertions(+), 128 deletions(-)

diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index 52b7a703f..d14474df0 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -34,8 +34,8 @@ void export_entity_to_density();
 
 namespace {
   
-Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real, const String&)=&mol::alg::LocalDistTest;
-Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
+Real (*ldt_a)(const mol::EntityView&, const mol::alg::GlobalDistanceList& glob_dist_list, Real, const String&)=&mol::alg::LocalDistTest;
+Real (*ldt_b)(const seq::AlignmentHandle&,Real, int, int)=&mol::alg::LocalDistTest;
 mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes;
 mol::EntityView (*fc_b)(const mol::EntityHandle&, const mol::alg::ClashingDistances&, bool)=&mol::alg::FilterClashes;
 mol::EntityView (*csc_a)(const mol::EntityView&, const mol::alg::StereoChemicalParams&, const mol::alg::StereoChemicalParams&, Real, Real, bool)=&mol::alg::CheckStereoChemistry;
diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc
index 6b343838c..3b37839c6 100644
--- a/modules/mol/alg/src/ldt.cc
+++ b/modules/mol/alg/src/ldt.cc
@@ -235,7 +235,8 @@ int main (int argc, char **argv)
       v=alg::FilterClashes(v,nonbonded_table);
     }
     
-    Real ldt=LDTHA(v, ref_view, radius);
+    GlobalDistanceList glob_dist_list = CreateDistanceList(ref_view,radius);
+    Real ldt=LDTHA(v, glob_dist_list);
     
     std::cout << "File: " << files[i] << std::endl; 
     std::cout << "Global LDT score: " << ldt << std::endl;
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
index 274a32aa1..13f37e2e0 100644
--- a/modules/mol/alg/src/local_dist_test.cc
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -1,6 +1,7 @@
 #include <ost/log.hh>
 #include <ost/mol/mol.hh>
 #include "local_dist_test.hh"
+#include <boost/concept_check.hpp>
 
 namespace ost { namespace mol { namespace alg {
 
@@ -37,7 +38,8 @@ bool swappable(const String& rname, const String& aname)
   if (rname=="ASP") {
     return (aname=="OD1" || aname=="OD2");
   }
-  if (rname=="VAL") {
+  if (rname=="VAL") {     
+
     return (aname=="CG1" || aname=="CG2");
   }
   if (rname=="TYR" || rname=="PHE") {
@@ -49,84 +51,103 @@ bool swappable(const String& rname, const String& aname)
   return false;
 }
 
-std::pair<Real, Real> calc_overlap1(ResidueView ref_res, 
-                                    EntityView ref,
+std::pair<bool,Real> within_tolerance(Real mdl_dist, const ReferenceDistance& ref_dist, Real tol)
+{
+  Real min = ref_dist.GetMaxDistance();
+  Real max = ref_dist.GetMaxDistance();
+  bool within_tol = false;
+  Real difference = 0;
+  if (mdl_dist>=min && mdl_dist <=max) {
+    within_tol=true;
+  } else if (mdl_dist < min && abs(min-mdl_dist) < tol) {
+    within_tol = true;
+  } else if (mdl_dist > max && abs(mdl_dist-max) < tol) {
+    within_tol = true;
+  }
+  if (within_tol == false) {
+    if (mdl_dist > max) {
+      difference = mdl_dist-(max+tol);
+    } else {
+      difference = mdl_dist-(min-tol);
+    }
+  }  
+  return std::make_pair<bool,Real>(within_tol,difference);  
+}      
+
+
+std::pair<Real, Real> calc_overlap1(const ResidueDistanceList& res_distance_list, 
                                     ChainView mdl_chain, 
-                                    Real tol, Real max_dist, bool only_fixed, 
-                                    bool swap,std::vector<std::pair<Real, Real> >& overlap_list )
+                                    Real tol, bool only_fixed, 
+                                    bool swap,std::vector<std::pair<Real, Real> >& overlap_list, bool log )
 {
   std::pair<Real, Real> overlap(0.0, 0.0);
-  AtomViewList ref_atoms=ref_res.GetAtomList();
 
-  ResidueView mdl_res=mdl_chain.FindResidue(ref_res.GetNumber());
-  AtomViewList within;
-  if (max_dist<0){
-      within=ref.GetAtomList();
-  }
-  for (AtomViewList::iterator ai=ref_atoms.begin(),
-       ae=ref_atoms.end(); ai!=ae; ++ai) {
-    if (ai->GetElement()=="H") { continue; }
-    String name=swap ? swapped_name(ai->GetName()) : ai->GetName();
+  ResidueView mdl_res=mdl_chain.FindResidue(res_distance_list[0].GetFirstAtom().GetResNum()); 
+  for (ResidueDistanceList::const_iterator ai=res_distance_list.begin(),
+       ae=res_distance_list.end(); ai!=ae; ++ai) {
+    UniqueAtomIdentifier first_atom=ai->GetFirstAtom();
+    UniqueAtomIdentifier second_atom=ai->GetSecondAtom();
+    String name=swap ? swapped_name(first_atom.GetAtomName()) : first_atom.GetAtomName();
     AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
-    if (max_dist>=0){ 
-      within=ref.FindWithin(ai->GetPos(), max_dist);
-    }
-    for (AtomViewList::iterator aj=within.begin(),
-         ae2=within.end(); aj!=ae2; ++aj) {
-      if (aj->GetElement()=="H" ||
-          aj->GetResidue().GetChain()!=ai->GetResidue().GetChain()) {
+ 
+    if (only_fixed) {
+      if (swappable(second_atom.GetResidueName(), second_atom.GetAtomName())) {
           continue;
       }
-
-      if (only_fixed) {
-        if (aj->GetResidue().GetNumber()==ref_res.GetNumber()) {
-          continue;
-        }
-        if (swappable(aj->GetResidue().GetName(), aj->GetName())) {
-          continue;
-        }
-        AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(), 
-                                        aj->GetName());
-        overlap.second+=1.0;
-        if (!(av1 && av2)) {
-          continue;
-        }
-        Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
-        Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
-        if (std::abs(mdl_dist-ref_dist)<tol) {
-          overlap.first+=1;
-        }
+      AtomView av2=mdl_chain.FindAtom(second_atom.GetResNum(),second_atom.GetAtomName());
+      overlap.second+=1.0;
+      if (!(av1 && av2)) {
         continue;
       }
-      if (aj->GetResidue().GetNumber()>ref_res.GetNumber()) {
-        AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(), 
-                                             aj->GetName());
-        overlap.second+=1.0;
-        overlap_list[ref_res.GetIndex()].second+=1.0;
-        overlap_list[aj->GetResidue().GetIndex()].second+=1.0;
-        if (!(av1 && av2)) {
-          continue;
-        }
-        Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
-        Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
-        if (std::abs(mdl_dist-ref_dist)<tol) {
+      Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+      ReferenceDistance ref_dist=*ai;
+      if (within_tolerance(mdl_dist,ref_dist,tol).first) {
+        if (log) {
           LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName() 
              << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " " 
-             << mdl_dist << " " << ref_dist << " " << mdl_dist-ref_dist << " " << tol << " " << "PASS")
-          overlap.first+=1;
-          overlap_list[ref_res.GetIndex()].first+=1.0;
-          overlap_list[aj->GetResidue().GetIndex()].first+=1.0;
-        } else {
+             << mdl_dist << " " << ref_dist.GetMinDistance() << " " << ref_dist.GetMaxDistance() << " " << tol << " " << "PASS")
+        }   
+        overlap.first+=1;
+      } else {
+        if (log) {
           LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName() 
              << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " " 
-             << mdl_dist << " " << ref_dist << " " << mdl_dist-ref_dist << " " << tol << " " << "FAIL")
-        } 
+             << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "FAIL")
+        }  
       }
-    }      
-  }
+      continue;
+    }
+    AtomView av2=mdl_chain.FindAtom(second_atom.GetResNum(), second_atom.GetAtomName());
+    overlap.second+=1.0;
+    if (av1) {
+      overlap_list[av1.GetResidue().GetIndex()].second+=1.0;
+    }
+    if (av2) {
+      overlap_list[av2.GetResidue().GetIndex()].second+=1.0;
+     }  
+    if (!(av1 && av2)) {
+      continue;
+    }
+    
+    Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+    ReferenceDistance ref_dist=*ai;
+    if (within_tolerance(mdl_dist,ref_dist,tol).first) {
+      LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName() 
+           << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " " 
+           << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "PASS")
+      overlap.first+=1;
+      overlap_list[av1.GetResidue().GetIndex()].first+=1.0;
+      overlap_list[av2.GetResidue().GetIndex()].first+=1.0;
+    } else {
+      LOG_VERBOSE("LDT:" << " " << av1.GetResidue().GetChain() << " " << av1.GetResidue().GetName() << " " << av1.GetResidue().GetNumber() << " " << av1.GetName() 
+             << " " << av2.GetResidue().GetChain() << " " << av2.GetResidue().GetName() << " " << av2.GetResidue().GetNumber() << " " << av2.GetName() << " " 
+             << mdl_dist << " " << ref_dist.GetMinDistance() << " " << mdl_dist-ref_dist.GetMinDistance() << " " << tol << " " << "FAIL")
+    } 
+  }    
   return overlap;
 }
 
+
 std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
                                     const seq::ConstSequenceHandle& mdl_seq,
                                     int pos, Real tol, Real max_dist, 
@@ -186,7 +207,8 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
         continue;
       } else {
         if (aj->GetResidue().GetNumber()>ref_res.GetNumber()) {
-          overlap.second+=1.0;
+          overlap.second+=1.0;    
+
           try {
            int aln_pos=ref_seq.GetPos(aj->GetResidue().GetIndex());
             ResidueView r2=mdl_seq.GetResidue(aln_pos);
@@ -212,81 +234,141 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
 
 }
 
-Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
-                   Real cutoff, Real max_dist, const String& local_ldt_property_string)
+bool ReferenceDistance::IsValid() const 
+{
+  if (mind_ == -1.0 and maxd_ == -1.0) {
+    return false;
+  } 
+  return true;
+}
+
+void ReferenceDistance::Print() const 
+{ 
+  if (this->IsValid() == true) {
+    std::cout << first_atom_.GetChainName() << " " << first_atom_.GetResNum() << " " << first_atom_.GetResidueName() << " " << first_atom_.GetAtomName() << " " <<
+                 second_atom_.GetChainName() << " " << second_atom_.GetResNum() << " " << second_atom_.GetResidueName() << " " << second_atom_.GetAtomName() << " " <<
+                 mind_ << " " << maxd_ << std::endl; 
+  } else {
+    std::cout << first_atom_.GetChainName() << " " << first_atom_.GetResNum() << " " << first_atom_.GetResidueName() << " Placeholder" << std::endl;
+  }
+}
+
+GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist)
+{
+ GlobalDistanceList dist_list; 
+ ResidueViewList ref_residues=ref.GetChainList()[0].GetResidueList();
+ for (ResidueViewList::iterator i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
+   ResidueView rview = (*i);
+//#   if (rview.IsPeptideLinking()) {
+   if (true) {
+     ResidueDistanceList res_dist_list;
+     AtomViewList ref_atoms=(*i).GetAtomList();
+     AtomViewList within;
+     if (max_dist<0){
+     dist_list.push_back(res_dist_list);
+       within=ref.GetAtomList();
+     }  
+     for (AtomViewList::iterator ai=ref_atoms.begin(), ae=ref_atoms.end(); ai!=ae; ++ai) {
+       UniqueAtomIdentifier first_atom(ai->GetResidue().GetChain().GetName(),ai->GetResidue().GetNumber(),ai->GetResidue().GetName(),ai->GetName());
+       if (ai->GetElement()=="H") { continue; }
+       if (max_dist>=0){ 
+         within=ref.FindWithin(ai->GetPos(), max_dist);
+       }
+       for (AtomViewList::iterator aj=within.begin(), ae2=within.end(); aj!=ae2; ++aj) {      
+         UniqueAtomIdentifier second_atom(aj->GetResidue().GetChain().GetName(),aj->GetResidue().GetNumber(),aj->GetResidue().GetName(),aj->GetName());
+         if (aj->GetElement()=="H" ||
+             aj->GetResidue().GetChain()!=ai->GetResidue().GetChain()) {
+             continue;
+         }
+         if (aj->GetResidue().GetNumber()>i->GetNumber()) {
+           Real dist=geom::Length(ai->GetPos()-aj->GetPos());
+           ReferenceDistance ref_dist(first_atom,second_atom,dist,dist);
+           res_dist_list.push_back(ref_dist);
+         }
+       }
+     }
+     if (res_dist_list.size()==0) {
+       UniqueAtomIdentifier current_residue_fake_atom(rview.GetChain().GetName(),rview.GetNumber(),rview.GetName(),"CA");
+       ReferenceDistance fake_ref_distance(current_residue_fake_atom,current_residue_fake_atom,-1.0,-1.0);
+       res_dist_list.push_back(fake_ref_distance);
+     }  
+     dist_list.push_back(res_dist_list);
+   }
+ }  
+ return dist_list;
+} 
+
+
+Real LocalDistTest(const EntityView& mdl, const GlobalDistanceList& glob_dist_list,
+                   Real cutoff, const String& local_ldt_property_string)
 {
-    
   if (!mdl.GetResidueCount()) {
     LOG_WARNING("model structures doesn't contain any residues");
     return 0.0;
   }
-
-  if (!ref.GetResidueCount()) {
-    LOG_WARNING("reference structures doesn't contain any residues");
+  if (glob_dist_list.size()==0) {
+    LOG_WARNING("global reference list is empty");
     return 0.0;
   }
-  ResidueViewList ref_residues=ref.GetChainList()[0].GetResidueList();
-  std::vector<std::pair<Real, Real> > overlap_list(ref_residues.size(), std::pair<Real, Real>(0.0, 0.0));
+  std::vector<std::pair<Real, Real> > overlap_list(mdl.GetResidueCount(), std::pair<Real, Real>(0.0, 0.0));
   ChainView mdl_chain=mdl.GetChainList()[0];  
   // Residues with symmetric side-chains require special treatment as there are 
   // two possible ways to name the atoms. We calculate the overlap with the 
   // fixed atoms and take the solution that gives bigger scores.
-  XCSEditor edi=ref.GetHandle().EditXCS(BUFFERED_EDIT);
-  for (ResidueViewList::iterator
-       i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
-    const String rname=i->GetName();
-    ResidueView mdl_res=mdl_chain.FindResidue(i->GetNumber());
-    if (!mdl_res) {
-     continue;
-    }     
-    if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || 
-         rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+  XCSEditor edi=mdl.GetHandle().EditXCS(BUFFERED_EDIT);
+  for (GlobalDistanceList::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
+    ResidueDistanceList rdl = *i;
+    if (!rdl[0].IsValid()) {
       continue;
     }
-
-    std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain, 
-                                          cutoff, max_dist, true, 
-                                          false, overlap_list);
-    std::pair<Real, Real> ov2=calc_overlap1(*i, ref, mdl_chain, 
-                                          cutoff, max_dist, true, 
-                                          true, overlap_list); 
+    
+    ResNum rnum =  rdl[0].GetFirstAtom().GetResNum() ;
+    String rname=rdl[0].GetFirstAtom().GetResidueName();
+    ResidueView mdl_res=mdl_chain.FindResidue(rnum);  
+    if (!mdl_res) {   
+      continue;
+    }
+    if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" || rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+      continue;
+    }
+    std::pair<Real, Real> ov1=calc_overlap1(*i, mdl_chain, 
+                                          cutoff, true, 
+                                          false, overlap_list,false);
+ 
+    std::pair<Real, Real> ov2=calc_overlap1(*i, mdl_chain, 
+                                          cutoff, true, 
+                                          true, overlap_list,false); 
     if (ov1.first/ov1.second<ov2.first/ov2.second) {
-     AtomViewList atoms=mdl_res.GetAtomList();
-     for (AtomViewList::iterator j=atoms.begin(), 
-          e2=atoms.end(); j!=e2; ++j) {
-       if (swappable(rname, j->GetName())) {
-         edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
-       }
-     }
+      AtomViewList atoms=mdl_res.GetAtomList();
+      for (AtomViewList::iterator j=atoms.begin(), e2=atoms.end(); j!=e2; ++j) {
+        if (swappable(rname, j->GetName())) {
+          edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
+        }
+      }
     }
   }
+  overlap_list.clear();
   std::pair<Real, Real> total_ov(0.0, 0.0);
-  for (ResidueViewList::iterator
-       i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
-     std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain, cutoff, 
-                                            max_dist, false, false, overlap_list);
-     total_ov.first+=ov1.first;
-     total_ov.second+=ov1.second;
-  }
-
-  // attach local per-residue LDT as property:
-  if(local_ldt_property_string!="") {
-    for (ResidueViewList::iterator
-       i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
-     std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain, cutoff, 
-                                            max_dist, false, false, overlap_list);
-     total_ov.first+=ov1.first;
-     total_ov.second+=ov1.second;
-     Real local_ldt=overlap_list[(*i).GetIndex()].first/(overlap_list[(*i).GetIndex()].second ? overlap_list[(*i).GetIndex()].second : 1);
-     if (mdl_chain.FindResidue((*i).GetNumber()).IsValid()) {
-      ResidueView mdl_res=mdl_chain.FindResidue((*i).GetNumber());
-      mdl_res.SetFloatProp(local_ldt_property_string, local_ldt);
-     }
+  for (GlobalDistanceList::const_iterator i=glob_dist_list.begin(), e=glob_dist_list.end(); i!=e; ++i) {
+    ResidueDistanceList rdl = *i;
+    if (rdl[0].IsValid()) {
+      std::pair<Real, Real> ov1=calc_overlap1(*i, mdl_chain, cutoff, 
+                                            false, false, overlap_list,true);
+      total_ov.first+=ov1.first;
+      total_ov.second+=ov1.second;       
     }
-      mdl_chain.SetFloatProp(local_ldt_property_string, total_ov.first/(total_ov.second ? total_ov.second : 1));
+    if(local_ldt_property_string!="") {
+      ResNum rn = rdl[0].GetFirstAtom().GetResNum();
+      ResidueView mdlr=mdl_chain.FindResidue(rn);  
+      if (mdlr.IsValid()) {
+        int mdl_res_index =mdlr.GetIndex();
+        Real local_ldt=overlap_list[mdl_res_index].first/(overlap_list[mdl_res_index].second ? overlap_list[mdl_res_index].second : 1);
+        ResidueView res_to_wr = mdl_chain.FindResidue((*i)[0].GetFirstAtom().GetResNum()); 
+        res_to_wr.SetFloatProp(local_ldt_property_string, local_ldt);
+      }  
+    }       
   }
   overlap_list.clear();
-
   return total_ov.first/(total_ov.second ? total_ov.second : 1);
 }
 
@@ -341,16 +423,17 @@ Real LocalDistTest(const ost::seq::AlignmentHandle& aln,
   return 0.0;
 }
 
-Real LDTHA(EntityView&v, const EntityView& ref_view, Real radius)
+Real LDTHA(EntityView& v, const GlobalDistanceList& global_dist_list)
 {
+    
     Real cutoffs[]={0.5,1,2,4};
     String labels[]={"localldt0.5","localldt1","localldt2","ldtlocal4"};
     Real ldt=0.0;   
     for (int n=0; n<4; ++n) { 
-      ldt+=alg::LocalDistTest(v, ref_view, cutoffs[n], radius,labels[n]);
+      ldt+=alg::LocalDistTest(v, global_dist_list, cutoffs[n], labels[n]);
     }      
     ldt/=4.0;    
     return ldt;
 }
 
-}}}
\ No newline at end of file
+}}}
diff --git a/modules/mol/alg/src/local_dist_test.hh b/modules/mol/alg/src/local_dist_test.hh
index c6dd4d7e0..74473ba29 100644
--- a/modules/mol/alg/src/local_dist_test.hh
+++ b/modules/mol/alg/src/local_dist_test.hh
@@ -25,17 +25,66 @@
 
 namespace ost { namespace mol { namespace alg {
   
+  
+class UniqueAtomIdentifier
+{
+  
+public:
+  UniqueAtomIdentifier(const String& chain,const ResNum& residue,const String& residue_name, const String& atom): chain_(chain),residue_(residue),residue_name_(residue_name),atom_(atom) {}  
+  
+  String GetChainName() const { return chain_; } 
+  ResNum GetResNum() const { return residue_; }  
+  String GetResidueName() const { return residue_name_; }
+  String GetAtomName() const { return atom_; }
+  
+private:
+
+  String chain_;
+  ResNum residue_;
+  String residue_name_;
+  String atom_;    
+};
+
+class ReferenceDistance
+{
+  
+public:
+  
+  ReferenceDistance(UniqueAtomIdentifier first_atom, UniqueAtomIdentifier second_atom, Real mind, Real maxd):
+       first_atom_(first_atom),second_atom_(second_atom),mind_(mind),maxd_(maxd) {}
+  UniqueAtomIdentifier GetFirstAtom() const {return first_atom_;}
+  UniqueAtomIdentifier GetSecondAtom() const {return second_atom_;}
+  Real GetMinDistance() const {return mind_ ;}
+  Real GetMaxDistance() const {return maxd_ ;}
+  bool IsValid() const; 
+  void Print() const; 
+  
+private:
+  
+  UniqueAtomIdentifier first_atom_;
+  UniqueAtomIdentifier second_atom_;
+  Real mind_;
+  Real maxd_;
+  bool valid;
+ 
+};
+
+typedef std::vector<ReferenceDistance> ResidueDistanceList;
+typedef std::vector<ResidueDistanceList> GlobalDistanceList;
+  
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
-                                         const EntityView& ref,
-                                         Real cutoff, Real max_dist,
+                                         const GlobalDistanceList& dist_list,
+                                         Real cutoff, 
                                          const String& local_ldt_property_string="");
 
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
-                                         Real cutoff, Real max_dist, 
+                                         Real cutoff, 
                                          int ref_index=0, int mdl_index=1);
 
 
-Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView&v, const EntityView& ref_view, Real radius);
+Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView& v, const GlobalDistanceList& global_dist_list);
+
+GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist);
 
 
 }}}
-- 
GitLab