added statistical functions

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2332 5a81b35b-ba03-0410-adc8-b2c5c5119f08

modules/base/pymod/CMakeLists.txt

@@ -7,4 +7,4 @@ set(OST_BASE_PYMOD_SOURCES
 
 pymod(NAME base OUTPUT_DIR ost 
       CPP ${OST_BASE_PYMOD_SOURCES} 
-      PY __init__.py settings.py)
+      PY __init__.py settings.py stutil.py)

modules/base/pymod/__init__.py

@@ -17,6 +17,7 @@
 # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 #------------------------------------------------------------------------------
 from _base import *
+from stutil import *
 
 from ost import geom
 from ost import io

modules/base/pymod/stutil.py

+import math
+from ost import mol
+
+def FloatValueExtract(func):
+  """
+  Decorator to wrap functions that take a list of float values. In addition to 
+  passing in a list of float values directly, it is possible to extract the 
+  values from attributes or generic properties.
+  """
+  def _dec(xs, prop=None, attr=None):
+    if prop!=None:
+      if len(xs)==0:
+        return func([])
+      assert attr==None
+      level=mol.Prop.Level.UNSPECIFIED
+      if isinstance(xs[0], mol.AtomBase):
+        level=mol.Prop.Level.ATOM
+      elif isinstance(xs[0], mol.ResidueBase):
+        level=mol.Prop.Level.RESIDUE
+      elif isinstance(xs[0], mol.ChainBase):
+        level=mol.Prop.Level.CHAIN
+      epm=mol.EntityPropertyMapper(prop, level)
+      vals=[]
+      for x in xs:
+        try:
+          vals.append(epm.Get(x))
+        except:
+          pass
+      return func(vals)
+    if attr!=None:
+      vals=[]
+      for x in xs:
+        try:
+          vals.append(getattr(x, attr))
+        except:
+          pass
+      return func(vals)
+    return func(xs)
+  return _dec
+
+@FloatValueExtract  
+def Mean(xs):
+  """
+  Calculate mean of dataset
+  """
+  if len(xs)==0:
+    raise RuntimeError("Can't calculate mean of empty sequence")
+  return sum(xs)/len(xs)
+
+@FloatValueExtract
+def StdDev(xs):
+  """
+  Calculate standard-deviation of dataset
+  
+            | sum[xi-<x>]^2 |
+  sigma=sqrt|---------------|
+            |       n       |
+  """
+  mean=Mean(xs)
+  return math.sqrt(sum([(x-mean)**2 for x in xs])/len(xs))
+
+@FloatValueExtract
+def Min(xs):
+  return min(xs)
+
+@FloatValueExtract
+def Max(xs):
+  return max(xs)
+
+def Correl(xs, ys):
+  """
+  Calculates the correlation coefficient between xs and ys as
+  
+    sum[(xi-<x>)*(yi-<y>)]
+  r=----------------------
+          (n-1)*sx*sy
+          
+  where <x>, <y> are the mean of dataset xs and ys, and, sx and sy are the 
+  standard deviations.
+  """
+  if len(xs)!=len(ys):
+    raise RuntimeError("Can't calculate correl. Sequence lengths do not match.")
+  if len(xs)==1:
+    raise RuntimeError("Can't calculate correl of sequences with length 1.")
+  mean_x=Mean(xs)
+  mean_y=Mean(ys)
+  sigma_x, sigma_y=(0.0, 0.0)
+  cross_term=0.0
+  for x, y in zip(xs, ys):
+    cross_term+=(x-mean_x)*(y-mean_y)
+    sigma_x+=(x-mean_x)**2
+    sigma_y+=(y-mean_y)**2
+  sigma_x=math.sqrt(sigma_x/len(xs))
+  sigma_y=math.sqrt(sigma_y/len(ys))  
+  return cross_term/((len(xs)-1)*sigma_x*sigma_y)

modules/mol/base/pymod/CMakeLists.txt

@@ -18,6 +18,7 @@ export_query_view_wrapper.cc
 export_torsion.cc
 export_visitor.cc
 wrap_mol.cc
+export_entity_property_mapper.cc
 )
 
 pymod(NAME mol CPP ${OST_BASE_PYMOD_SOURCES} PY __init__.py)

modules/mol/base/pymod/export_entity_property_mapper.cc

+#include <boost/python.hpp>
+
+#include <ost/mol/entity_property_mapper.hh>
+
+using namespace ost;
+using namespace ost::mol;
+using namespace boost::python;
+
+namespace {
+
+
+Real (EntityPropertyMapper::*get_ah_a)(const AtomHandle&) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_av_a)(const AtomView&) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_rh_a)(const ResidueHandle&) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_rv_a)(const ResidueView&) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_ch_a)(const ChainHandle&) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_cv_a)(const ChainView&) const=&EntityPropertyMapper::Get;
+
+
+Real (EntityPropertyMapper::*get_ah_b)(const AtomHandle&, Real) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_av_b)(const AtomView&, Real) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_rh_b)(const ResidueHandle&, Real) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_rv_b)(const ResidueView&, Real) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_ch_b)(const ChainHandle&, Real) const=&EntityPropertyMapper::Get;
+Real (EntityPropertyMapper::*get_cv_b)(const ChainView&, Real) const=&EntityPropertyMapper::Get;
+
+EntityPropertyMapper create_epm(const String& prop_name, char level)
+{
+  switch(level) {
+    case 'A':
+    case 'a':
+      return EntityPropertyMapper(prop_name, Prop::ATOM);
+    case 'R':
+    case 'r':
+      return EntityPropertyMapper(prop_name, Prop::RESIDUE);
+    case 'C':
+    case 'c':
+      return EntityPropertyMapper(prop_name, Prop::CHAIN);
+    case 'U':
+    case 'u':
+      return EntityPropertyMapper(prop_name, Prop::UNSPECIFIED);
+    default:
+      throw std::runtime_error(String("unknown property level '")+level+"'");
+  }
+}
+
+}
+
+void export_EntityPropertyMapper()
+{
+  class_<EntityPropertyMapper>("EntityPropertyMapper", 
+                               init<const String&, 
+                                    Prop::Level>(arg("level")=Prop::UNSPECIFIED))
+    .def("__init__", &create_epm)
+    .def("Get", get_ah_a)
+    .def("Get", get_av_a)
+    .def("Get", get_rh_a)
+    .def("Get", get_rv_a)
+    .def("Get", get_ch_a)
+    .def("Get", get_cv_a)
+    .def("Get", get_ah_b)
+    .def("Get", get_av_b)
+    .def("Get", get_rh_b)
+    .def("Get", get_rv_b)
+    .def("Get", get_ch_b)
+    .def("Get", get_cv_b)
+  ;  
+}

modules/mol/base/pymod/wrap_mol.cc

@@ -41,6 +41,7 @@ void export_CoordGroup();
 void export_PropertyID();
 void export_BoundingBox();
 void export_QueryViewWrapper();
+void export_EntityPropertyMapper();
 BOOST_PYTHON_MODULE(_mol)
 {
   export_Entity();
@@ -61,7 +62,7 @@ BOOST_PYTHON_MODULE(_mol)
   export_PropertyID();  
   export_BoundingBox();
   export_QueryViewWrapper();
-
+  export_EntityPropertyMapper();
   class_<Transform>("Transform", init<>())
     .def("GetMatrix",&Transform::GetMatrix)
     .def("GetTransposedMatrix",&Transform::GetTransposedMatrix)