Local per reside LDT

Adds to the LDT function the option to save the local LDT in a floatProp. Mainly work by Pascal

modules/mol/alg/pymod/wrap_mol_alg.cc

@@ -22,19 +22,18 @@
 #include <ost/mol/alg/local_dist_test.hh>
 #include <ost/mol/alg/superpose_frames.hh>
 #include <ost/mol/alg/filter_clashes.hh>
-#include <ost/mol/alg/construct_cbeta.hh>
 using namespace boost::python;
 using namespace ost;
 
 void export_svdSuperPose();
-void export_Clash();
+
 #if OST_IMG_ENABLED
 void export_entity_to_density();
 #endif
 
 namespace {
   
-Real (*ldt_a)(const mol::EntityView&, const mol::EntityView& ref, Real, Real)=&mol::alg::LocalDistTest;
+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;
 mol::EntityView (*fc_a)(const mol::EntityView&, Real,bool)=&mol::alg::FilterClashes;
 mol::EntityView (*fc_b)(const mol::EntityHandle&, Real, bool)=&mol::alg::FilterClashes;
@@ -47,15 +46,11 @@ BOOST_PYTHON_MODULE(_mol_alg)
   export_entity_to_density();
   #endif
   
-  def("LocalDistTest", ldt_a);
+  def("LocalDistTest", ldt_a, (arg("local_ldt_property_string")=""));
   def("LocalDistTest", ldt_b, (arg("ref_index")=0, arg("mdl_index")=1));
   def("FilterClashes", fc_a, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
   def("FilterClashes", fc_b, (arg("ent"), arg("tolerance")=0.1, arg("always_remove_bb")=false));
   def("SuperposeFrames", &ost::mol::alg::SuperposeFrames, 
       (arg("source"), arg("sel")=ost::mol::EntityView(), arg("begin")=0, 
        arg("end")=-1, arg("ref")=-1));
-//   def("ConstructCBetas", one_arg, args("entity_handle"));
-  def("ConstructCBetas", &ost::mol::alg::ConstructCBetas, (arg("entity_handle"), arg("include_gly")=false));
-  
-  export_Clash();
 }

modules/mol/alg/src/local_dist_test.cc

@@ -53,7 +53,7 @@ std::pair<Real, Real> calc_overlap1(ResidueView ref_res,
                                     EntityView ref,
                                     ChainView mdl_chain, 
                                     Real tol, Real max_dist, bool only_fixed, 
-                                    bool swap)
+                                    bool swap,std::vector<std::pair<Real, Real> >& overlap_list )
 {
   std::pair<Real, Real> overlap(0.0, 0.0);
   AtomViewList ref_atoms=ref_res.GetAtomList();
@@ -98,6 +98,8 @@ std::pair<Real, Real> calc_overlap1(ResidueView ref_res,
         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;
         }
@@ -105,6 +107,8 @@ std::pair<Real, Real> calc_overlap1(ResidueView ref_res,
         Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
         if (std::abs(mdl_dist-ref_dist)<tol) {
           overlap.first+=1;
+          overlap_list[ref_res.GetIndex()].first+=1.0;
+          overlap_list[aj->GetResidue().GetIndex()].first+=1.0;
         }
       }
     }      
@@ -195,8 +199,9 @@ std::pair<Real, Real> calc_overlap2(const seq::ConstSequenceHandle& ref_seq,
 }
 
 Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
-                   Real cutoff, Real max_dist)
+                   Real cutoff, Real max_dist, const String& local_ldt_property_string)
 {
+    
   if (!mdl.GetResidueCount()) {
     LOG_WARNING("model structures doesn't contain any residues");
     return 0.0;
@@ -207,6 +212,7 @@ Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
     return 0.0;
   }
   ResidueViewList ref_residues=ref.GetResidueList();  
+  std::vector<std::pair<Real, Real> > overlap_list(ref_residues.size(), 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 
@@ -226,10 +232,10 @@ Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
 
     std::pair<Real, Real> ov1=calc_overlap1(*i, ref, mdl_chain, 
                                           cutoff, max_dist, true, 
-                                          false);
+                                          false, overlap_list);
     std::pair<Real, Real> ov2=calc_overlap1(*i, ref, mdl_chain, 
                                           cutoff, max_dist, true, 
-                                          true); 
+                                          true, overlap_list); 
     if (ov1.first/ov1.second<ov2.first/ov2.second) {
      AtomViewList atoms=mdl_res.GetAtomList();
      for (AtomViewList::iterator j=atoms.begin(), 
@@ -244,10 +250,29 @@ Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
   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);
+                                            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);
+     }
+    }
+      mdl_chain.SetFloatProp(local_ldt_property_string, total_ov.first/(total_ov.second ? total_ov.second : 1));
+  }
+  overlap_list.clear();
+
   return total_ov.first/(total_ov.second ? total_ov.second : 1);
 }
 

modules/mol/alg/src/local_dist_test.hh

@@ -27,12 +27,15 @@ namespace ost { namespace mol { namespace alg {
   
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
                                          const EntityView& ref,
-                                         Real cutoff, Real max_dist);
+                                         Real cutoff, Real max_dist,
+                                         const String& local_ldt_property_string="");
 
 Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const ost::seq::AlignmentHandle& aln,
                                          Real cutoff, Real max_dist, 
                                          int ref_index=0, int mdl_index=1);
+
 }}}
 
 #endif
 
+