New map-based ldt implementation

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -18,7 +18,7 @@
 //------------------------------------------------------------------------------
 
 #include <boost/python.hpp>
-#include <boost/python/suite/indexing/vector_indexing_suite.hpp>
+#include <boost/python/suite/indexing/map_indexing_suite.hpp>
 #include <ost/config.hh>
 #include <ost/mol/alg/local_dist_test.hh>
 #include <ost/mol/alg/superpose_frames.hh>
@@ -36,7 +36,7 @@ void export_entity_to_density();
 
 namespace {
   
-Real (*ldt_a)(const mol::EntityView&, const mol::alg::GlobalDistanceList& , Real, const String&)=&mol::alg::LocalDistTest;
+Real (*ldt_a)(const mol::EntityView&, const mol::alg::GlobalRDMap& , std::vector<Real>, const String&)=&mol::alg::LocalDistTest;
 Real (*ldt_c)(const mol::EntityView&, const mol::EntityView& , Real, Real, const String&)=&mol::alg::LocalDistTest;
 Real (*ldt_b)(const seq::AlignmentHandle&,Real, Real, int, int)=&mol::alg::LocalDistTest;
 mol::EntityView (*fc_a)(const mol::EntityView&, const mol::alg::ClashingDistances&,bool)=&mol::alg::FilterClashes;
@@ -70,8 +70,19 @@ ost::mol::alg::ClashingDistances fill_clashing_distances_wrapper (const list& st
  return ost::mol::alg::FillClashingDistances(stereo_chemical_props_file_vector,min_default_distance,min_distance_tolerance);
 }
 
+ost::mol::alg::GlobalRDMap create_distance_list_from_multiple_references(const list& ref_list, Real cutoff, Real max_dist)
+{
+  int ref_list_length = boost::python::extract<int>(ref_list.attr("__len__")());
+  std::vector<ost::mol::EntityView> ref_list_vector(ref_list_length);
+  
+  for (int i=0; i<ref_list_length; i++) {
+    ref_list_vector[i] = boost::python::extract<ost::mol::EntityView>(ref_list[i]);
+  }
+ 
+  return ost::mol::alg::CreateDistanceListFromMultipleReferences(ref_list_vector, cutoff, max_dist);  	
 }
 
+}
 
 
 BOOST_PYTHON_MODULE(_ost_mol_alg)
@@ -92,6 +103,7 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
   def("CheckStereoChemistry", csc_b, (arg("ent"), arg("bonds"), arg("angles"), arg("bond_tolerance"), arg("angle_tolerance"), arg("always_remove_bb")=false));
   def("LDTHA",&mol::alg::LDTHA);
   def("CreateDistanceList",&mol::alg::CreateDistanceList);
+  def("CreateDistanceListFromMultipleReferences",&create_distance_list_from_multiple_references);
     
   def("SuperposeFrames", superpose_frames1, 
       (arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0, 
@@ -127,23 +139,18 @@ BOOST_PYTHON_MODULE(_ost_mol_alg)
   ;    
    
   
-  class_<mol::alg::ReferenceDistance> ("ReferenceDistance", init <const mol::alg::UniqueAtomIdentifier&,const mol::alg::UniqueAtomIdentifier&, Real, Real>())
-    .def("GetFirstAtom",&mol::alg::ReferenceDistance::GetFirstAtom)
-    .def("GetSecondAtom",&mol::alg::ReferenceDistance::GetSecondAtom)
-    .def("GetMinDistance",&mol::alg::ReferenceDistance::GetMinDistance)
-    .def("GetMaxDistance",&mol::alg::ReferenceDistance::GetMaxDistance)
-  ;
-  
-  class_<std::vector<mol::alg::ReferenceDistance> >("ResidueDistanceList")
-    .def(vector_indexing_suite<std::vector<mol::alg::ReferenceDistance > >())
+  class_<mol::alg::ResidueRDMap>("ResidueRDMap")
+    .def(map_indexing_suite<mol::alg::ResidueRDMap>())
   ;
   
-  class_<std::vector<mol::alg::ResidueDistanceList> >("GlobalDistanceList")
-    .def(vector_indexing_suite<std::vector<mol::alg::ResidueDistanceList > >())
+  class_<mol::alg::GlobalRDMap>("GlobalRDMap")
+    .def(map_indexing_suite<mol::alg::GlobalRDMap>())
   ;
   
   def("FillClashingDistances",&fill_clashing_distances_wrapper);
   def("FillStereoChemicalParams",&fill_stereochemical_params_wrapper);
   def("IsStandardResidue",&mol::alg::IsStandardResidue);
+  def("PrintGlobalRDMap",&mol::alg::PrintGlobalRDMap);
+  def("PrintResidueRDMap",&mol::alg::PrintResidueRDMap);
   
 }

modules/mol/alg/src/ldt.cc

@@ -77,7 +77,7 @@ void usage()
   std::cerr << "   -e         print version" << std::endl;
 }
 
-std::pair<int,int> compute_coverage (const EntityView& v,const GlobalDistanceList& glob_dist_list)
+std::pair<int,int> compute_coverage (const EntityView& v,const GlobalRDMap& glob_dist_list)
 {
   int second=0;
   int first=0;
@@ -85,15 +85,12 @@ std::pair<int,int> compute_coverage (const EntityView& v,const GlobalDistanceLis
     return std::make_pair<int,int>(0,1);
   }
   ChainView vchain=v.GetChainList()[0];
-  for (std::vector<ResidueDistanceList>::const_iterator i=glob_dist_list.begin();i!=glob_dist_list.end();++i)
+  for (GlobalRDMap::const_iterator i=glob_dist_list.begin();i!=glob_dist_list.end();++i)
   {
-    ResNum rnum = (*i)[0].GetFirstAtom().GetResNum();
-    String rname = (*i)[0].GetFirstAtom().GetResidueName();
-    if (IsStandardResidue(rname)) {
-      second++;
-      if (vchain.FindResidue(rnum)) {
-        first++;
-      }
+    ResNum rnum = (*i).first;
+    second++;
+    if (vchain.FindResidue(rnum)) {
+      first++;
     }
   }
   return std::make_pair<int,int>(first,second);
@@ -200,7 +197,7 @@ int main (int argc, char **argv)
   }
   files.pop_back();
   EntityView ref_view=ref.Select(sel);
-  GlobalDistanceList glob_dist_list = CreateDistanceList(ref_view,radius);
+  GlobalRDMap glob_dist_list = CreateDistanceList(ref_view,radius);
   std::cout << "Verbosity level: " << verbosity_level << std::endl;
   if (structural_checks) {
     std::cout << "Stereo-chemical and steric clash checks: On " << std::endl;
@@ -296,22 +293,18 @@ int main (int argc, char **argv)
       std::cout << "Global LDT score: 0.0" << std::endl;
       return 0;
     }
-    Real ldt=LDTHA(v, glob_dist_list);
+    Real ldt=LDTHA(v, glob_dist_list);	
 
     std::cout << "Global LDT score: " << ldt << std::endl;
     std::cout << "Local LDT Score:" << std::endl;
     std::cout << "Chain\tResName\tResNum\tScore" << std::endl;
-    String labels[]={"localldt0.5","localldt1","localldt2","ldtlocal4"};
     for (ResidueViewIter rit=v.ResiduesBegin();rit!=v.ResiduesEnd();++rit){
       ResidueView ritv = *rit;
-      Real ldt_local_sum = 0;
-      if (ritv.HasProp("localldt0.5")) {
-        for (int n=0; n<4; ++n) {
-          ldt_local_sum+=ritv.GetFloatProp(labels[n]);
-        }
-        ldt_local_sum/=4.0;
-        std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << ldt_local_sum << std::endl;
+      Real ldt_local = 0;
+      if (ritv.HasProp("localldt")) {
+          ldt_local=ritv.GetFloatProp("localldt");
       }
+      std::cout << ritv.GetChain() << "\t" << ritv.GetName() << "\t" << ritv.GetNumber() << '\t' << ldt_local << std::endl;
     }
     std::cout << std::endl;
   }

modules/mol/alg/src/local_dist_test.hh

@@ -31,12 +31,16 @@ 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) {}  
-  
+
+  // to make the compiler happy (boost python map suite)
+  UniqueAtomIdentifier(): chain_(""),residue_(ResNum(1)),residue_name_(""),atom_("") {}  
+    
   String GetChainName() const { return chain_; } 
   ResNum GetResNum() const { return residue_; }  
   String GetResidueName() const { return residue_name_; }
   String GetAtomName() const { return atom_; }
   bool operator==(const UniqueAtomIdentifier& rhs) const;
+  bool operator<(const UniqueAtomIdentifier& rhs) const;	
     
 private:
 
@@ -46,40 +50,19 @@ private:
   String atom_;    
 };
 
-class ReferenceDistance
-{
-  
-public:
-  
-  ReferenceDistance(const UniqueAtomIdentifier& first_atom, const 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; 
-  bool operator==(const ReferenceDistance& rhs) const;
-  
-private:
-  
-  UniqueAtomIdentifier first_atom_;
-  UniqueAtomIdentifier second_atom_;
-  Real mind_;
-  Real maxd_;
-};
+typedef std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier> UAtomIdentifiers;
+typedef std::map<std::pair<UniqueAtomIdentifier,UniqueAtomIdentifier>,std::pair<float,float> > ResidueRDMap;
+typedef std::map<ost::mol::ResNum,ResidueRDMap> GlobalRDMap;
+typedef std::map<UniqueAtomIdentifier,int> ExistenceMap;
 
-typedef std::vector<ReferenceDistance> ResidueDistanceList;
-typedef std::vector<ResidueDistanceList> GlobalDistanceList;
-  
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
-                                         const GlobalDistanceList& dist_list,
-                                         Real cutoff, 
+                                         const GlobalRDMap& dist_list,
+                                         std::vector<Real> cutoff_list, 
                                          const String& local_ldt_property_string="");
 
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
                                          const EntityView& target,
-                                         Real cutoff, 
+                                         Real cutoff_list, 
                                          Real max_dist,
                                          const String& local_ldt_property_string="");
 
@@ -88,12 +71,15 @@ Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
                                          int ref_index=0, int mdl_index=1);
 
 
-Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView& v, const GlobalDistanceList& global_dist_list);
-
-GlobalDistanceList CreateDistanceList(const EntityView& ref,Real max_dist);
-bool IsStandardResidue(String rn);
+Real DLLEXPORT_OST_MOL_ALG LDTHA(EntityView& v, const GlobalRDMap& global_dist_list);
+Real DLLEXPORT_OST_MOL_ALG OldStyleLDTHA(EntityView& v, const GlobalRDMap& global_dist_list);
 
 
+GlobalRDMap CreateDistanceList(const EntityView& ref,Real max_dist);
+GlobalRDMap CreateDistanceListFromMultipleReferences(const std::vector<EntityView>& ref_list,Real cutoff, Real max_dist);
+void PrintGlobalRDMap(const GlobalRDMap& glob_dist_list);
+void PrintResidueRDMap(const ResidueRDMap& res_dist_list);
+bool IsStandardResidue(String rn);
 
 }}}
 

modules/mol/base/src/residue_prop.hh

@@ -34,6 +34,12 @@ class DLLEXPORT ResNum: private
     boost::unit_steppable<ResNum> > > > >
 {
 public:
+ 
+  // needed to wrap certain map classes
+  ResNum():
+    num_(1),alt_('\0')	 
+  {}  
+
   ResNum(int n):
     num_(n), alt_('\0')
   { }