From 9d02e0991a4364ed09a3c2ba623c74378dde1e1e Mon Sep 17 00:00:00 2001 From: Ansgar Philippsen <ansgar.philippsen@gmail.com> Date: Mon, 21 Feb 2011 17:32:00 -0500 Subject: [PATCH] moved and expanded numpy setpos to xcs editor various related cleanups and optimizations: - removed superfluous sort in load_dcd - added explicit *TransformedPos() methods to complement *Pos() methods (XCSEditor, AtomHandle, AtomImpl) - changed getter/setter logic for AtomImpl Name, TransformedPos and OriginalPos to direct access, deprecated GetName and removed Get*Pos and Set*Pos --- modules/io/src/mol/dcd_io.cc | 6 +- modules/mol/base/pymod/export_editors.cc | 116 +++++++++++++++++++- modules/mol/base/pymod/export_entity.cc | 100 +++++------------ modules/mol/base/src/atom_base.cc | 8 +- modules/mol/base/src/editor_base.cc | 2 +- modules/mol/base/src/impl/atom_impl.cc | 12 +- modules/mol/base/src/impl/atom_impl.hh | 20 ++-- modules/mol/base/src/impl/chain_impl.cc | 8 +- modules/mol/base/src/impl/connector_impl.cc | 8 +- modules/mol/base/src/impl/dihedral.cc | 4 +- modules/mol/base/src/impl/entity_impl.cc | 31 +++--- modules/mol/base/src/impl/residue_impl.cc | 49 +++++---- modules/mol/base/src/query_state.cc | 10 +- modules/mol/base/src/xcs_editor.cc | 106 ++++++++++++++++-- modules/mol/base/src/xcs_editor.hh | 39 ++++++- modules/mol/base/tests/test_numpy.py | 17 ++- 16 files changed, 364 insertions(+), 172 deletions(-) diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc index 3fdb02f55..e5f48a93d 100644 --- a/modules/io/src/mol/dcd_io.cc +++ b/modules/io/src/mol/dcd_io.cc @@ -184,7 +184,7 @@ bool read_frame(std::istream& istream, const DCDHeader& header, } -mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2, +mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted! const String& trj_fn, unsigned int stride) { @@ -196,10 +196,6 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2, } Profile profile_load("LoadCHARMMTraj"); - mol::AtomHandleList alist(alist2); - std::sort(alist.begin(),alist.end(),less_index); - - DCDHeader header; bool swap_flag=false, skip_flag=false, gap_flag=false; read_dcd_header(istream, header, swap_flag, skip_flag, gap_flag); diff --git a/modules/mol/base/pymod/export_editors.cc b/modules/mol/base/pymod/export_editors.cc index 9976b0d8f..87c19665b 100644 --- a/modules/mol/base/pymod/export_editors.cc +++ b/modules/mol/base/pymod/export_editors.cc @@ -26,6 +26,10 @@ using namespace boost::python; using namespace ost; using namespace ost::mol; +#if OST_NUMPY_SUPPORT_ENABLED +#include <numpy/arrayobject.h> +#endif + namespace { BondHandle (EditorBase::*connect_a)(const AtomHandle&, @@ -53,11 +57,111 @@ void (ICSEditor::*rotate_torsion_a)(TorsionHandle, Real)=&ICSEditor::RotateTorsi void (ICSEditor::*rotate_torsion_b)(const AtomHandle&, const AtomHandle&, const AtomHandle&, const AtomHandle&, Real)=&ICSEditor::RotateTorsionAngle; - + +void (XCSEditor::*set_pos1)(const AtomHandle&, const geom::Vec3&) = &XCSEditor::SetAtomPos; +void (XCSEditor::*set_t_pos1)(const AtomHandle&, const geom::Vec3&) = &XCSEditor::SetAtomTransformedPos; +void (XCSEditor::*set_o_pos1)(const AtomHandle&, const geom::Vec3&) = &XCSEditor::SetAtomOriginalPos; + +#if OST_NUMPY_SUPPORT_ENABLED +template<typename T, bool O> +void set_pos2_nc_t(XCSEditor& e, const AtomHandleList& alist, PyArrayObject* na) +{ + size_t count=0; + for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait,++count) { + if(O) { + e.SetAtomOriginalPos(*ait,geom::Vec3(static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,0))), + static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,1))), + static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,2))))); + } else { + e.SetAtomTransformedPos(*ait,geom::Vec3(static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,0))), + static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,1))), + static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,2))))); + } + } +} + +template<bool O> +void set_pos2_t(XCSEditor& e, const AtomHandleList& alist, object pyobj) +{ + size_t acount = alist.size(); + + if(!PyArray_Check(pyobj.ptr())) { + throw std::runtime_error("expected a numpy array"); + return; + } + PyArrayObject* na=reinterpret_cast<PyArrayObject*>(pyobj.ptr()); + + if(PyArray_NDIM(na)!=2 || PyArray_DIM(na,0)!=acount || PyArray_DIM(na,1)!=3) { + throw std::runtime_error("excpected a numpy array of shape (NAtoms, 3)"); + return; + } + + if(PyArray_ISCONTIGUOUS(na)) { + if(PyArray_TYPE(na)==NPY_FLOAT) { + if(O) { + e.SetAtomOriginalPos(alist,reinterpret_cast<float*>(PyArray_DATA(na))); + } else { + e.SetAtomTransformedPos(alist,reinterpret_cast<float*>(PyArray_DATA(na))); + } + } else if(PyArray_TYPE(na)==NPY_DOUBLE) { + if(O) { + e.SetAtomOriginalPos(alist,reinterpret_cast<double*>(PyArray_DATA(na))); + } else { + e.SetAtomTransformedPos(alist,reinterpret_cast<double*>(PyArray_DATA(na))); + } + } else { + throw std::runtime_error("excpected a numpy array of type float or double"); + return; + } + } else { + // non-contiguous +#if 0 + throw std::runtime_error("excpected contiguous numpy array"); +#else + if(PyArray_TYPE(na)==NPY_FLOAT) { + set_pos2_nc_t<float,O>(e,alist,na); + } else if(PyArray_TYPE(na)==NPY_DOUBLE) { + set_pos2_nc_t<double,O>(e,alist,na); + } else { + throw std::runtime_error("excpected a numpy array of type float or double"); + return; + } +#endif + } +} +#endif + +void set_t_pos2(XCSEditor& e, const AtomHandleList& alist, object pyobj) +{ +#if OST_NUMPY_SUPPORT_ENABLED + set_pos2_t<false>(e,alist,pyobj); +#else + throw std::runtime_error("SetAtomTransformedPos(alist,ndarray) disabled, since numpy support is not compiled in"); +#endif +} + +void set_o_pos2(XCSEditor& e, const AtomHandleList& alist, object pyobj) +{ +#if OST_NUMPY_SUPPORT_ENABLED + set_pos2_t<true>(e,alist,pyobj); +#else + throw std::runtime_error("SetAtomOriginalPos(alist,ndarray) disabled, since numpy support is not compiled in"); +#endif +} + +void set_pos2(XCSEditor& e, const AtomHandleList& alist, object pyobj) +{ + set_t_pos2(e,alist,pyobj); +} + } void export_Editors() -{ +{ +#if OST_NUMPY_SUPPORT_ENABLED + import_array(); +#endif + class_<EditorBase>("EditorBase", no_init) .def("InsertChain", &EditorBase::InsertChain) .def("InsertAtom", &EditorBase::InsertAtom, @@ -82,8 +186,12 @@ void export_Editors() ; class_<XCSEditor, bases<EditorBase> >("XCSEditor", no_init) - .def("SetAtomPos", &XCSEditor::SetAtomPos) - .def("SetAtomOriginalPos", &XCSEditor::SetAtomOriginalPos) + .def("SetAtomPos", set_pos1) + .def("SetAtomPos", set_pos2) + .def("SetAtomTransformedPos", set_t_pos1) + .def("SetAtomTransformedPos", set_t_pos2) + .def("SetAtomOriginalPos", set_o_pos1) + .def("SetAtomOriginalPos", set_o_pos2) .def("ApplyTransform", &XCSEditor::ApplyTransform) .def("SetTransform", &XCSEditor::SetTransform) .def("UpdateICS", &XCSEditor::UpdateICS) diff --git a/modules/mol/base/pymod/export_entity.cc b/modules/mol/base/pymod/export_entity.cc index 36ab544f1..826a07bb0 100644 --- a/modules/mol/base/pymod/export_entity.cc +++ b/modules/mol/base/pymod/export_entity.cc @@ -68,18 +68,34 @@ ICSEditor depr_request_ics_editor(EntityHandle e, EditMode m) return e.EditICS(m); } +bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2) +{ + return a1.GetIndex()<a2.GetIndex(); +} -PyObject* get_pos(EntityHandle& entity) + +PyObject* get_pos2(EntityHandle& entity, bool id_sorted) { #if OST_NUMPY_SUPPORT_ENABLED npy_intp dims[]={entity.GetAtomCount(),3}; - PyObject* na = PyArray_SimpleNew(2,dims,NPY_DOUBLE); - npy_double* nad = reinterpret_cast<npy_double*>(PyArray_DATA(na)); - for(AtomHandleIter it=entity.AtomsBegin();it!=entity.AtomsEnd();++it,nad+=3) { - geom::Vec3 pos=(*it).GetPos(); - nad[0]=static_cast<npy_double>(pos[0]); - nad[1]=static_cast<npy_double>(pos[1]); - nad[2]=static_cast<npy_double>(pos[2]); + PyObject* na = PyArray_SimpleNew(2,dims,NPY_FLOAT); + npy_float* nad = reinterpret_cast<npy_float*>(PyArray_DATA(na)); + if(id_sorted) { + AtomHandleList alist = entity.GetAtomList(); + std::sort(alist.begin(),alist.end(),less_index); + for(AtomHandleList::const_iterator it=alist.begin();it!=alist.end();++it,nad+=3) { + geom::Vec3 pos=(*it).GetPos(); + nad[0]=static_cast<npy_float>(pos[0]); + nad[1]=static_cast<npy_float>(pos[1]); + nad[2]=static_cast<npy_float>(pos[2]); + } + } else { + for(AtomHandleIter it=entity.AtomsBegin();it!=entity.AtomsEnd();++it,nad+=3) { + geom::Vec3 pos=(*it).GetPos(); + nad[0]=static_cast<npy_float>(pos[0]); + nad[1]=static_cast<npy_float>(pos[1]); + nad[2]=static_cast<npy_float>(pos[2]); + } } return na; #else @@ -88,67 +104,9 @@ PyObject* get_pos(EntityHandle& entity) #endif } - -#if OST_NUMPY_SUPPORT_ENABLED -template <typename T> -void set_pos_t(PyArrayObject* na, EntityHandle& entity) +PyObject* get_pos1(EntityHandle& entity) { - XCSEditor ed=entity.EditXCS(BUFFERED_EDIT); - - if(PyArray_ISCONTIGUOUS(na)) { - T* data = reinterpret_cast<T*>(PyArray_DATA(na)); - size_t count=0; - for(AtomHandleIter it=entity.AtomsBegin();it!=entity.AtomsEnd();++it,++count) { - ed.SetAtomPos(*it,geom::Vec3(static_cast<Real>(data[count*3+0]), - static_cast<Real>(data[count*3+1]), - static_cast<Real>(data[count*3+2]))); - } - } else { - size_t count=0; - for(AtomHandleIter it=entity.AtomsBegin();it!=entity.AtomsEnd();++it,++count) { - ed.SetAtomPos(*it,geom::Vec3(static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,0))), - static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,1))), - static_cast<Real>(*reinterpret_cast<T*>(PyArray_GETPTR2(na,count,2))))); - } - } -} -#endif - -void set_pos(EntityHandle& entity, object& pyobj) -{ -#if OST_NUMPY_SUPPORT_ENABLED - size_t acount = entity.GetAtomCount(); - - if(!PyArray_Check(pyobj.ptr())) { - throw std::runtime_error(std::string("expected a numpy array")); - return; - } - PyArrayObject* na=reinterpret_cast<PyArrayObject*>(pyobj.ptr()); - - if(PyArray_NDIM(na)!=2 || PyArray_DIM(na,0)!=acount || PyArray_DIM(na,1)!=3) { - throw std::runtime_error("excpected a numpy array of shape (NAtoms, 3)"); - return; - } - - switch(PyArray_TYPE(na)) { - case NPY_FLOAT: - set_pos_t<float>(na,entity); break; - case NPY_DOUBLE: - set_pos_t<double>(na,entity); break; - case NPY_INT: - set_pos_t<int>(na,entity); break; - case NPY_LONG: - set_pos_t<long>(na,entity); break; - case NPY_SHORT: - set_pos_t<short>(na,entity); break; - default: - throw std::runtime_error("excpected a numpy array of type float, double, int, long or short"); - return; - }; - -#else - throw std::runtime_error("SetPositions disabled, since numpy support is not compiled in"); -#endif + return get_pos2(entity,true); } } // ns @@ -230,9 +188,9 @@ void export_Entity() .def(self==self) .def(self!=self) #if OST_NUMPY_SUPPORT_ENABLED - .def("SetPositions",set_pos) - .def("GetPositions",get_pos) - .add_property("positions",get_pos,set_pos) + .def("GetPositions",get_pos1) + .def("GetPositions",get_pos2) + .add_property("positions",get_pos1) #endif ; diff --git a/modules/mol/base/src/atom_base.cc b/modules/mol/base/src/atom_base.cc index f25d40135..642ea5aa9 100644 --- a/modules/mol/base/src/atom_base.cc +++ b/modules/mol/base/src/atom_base.cc @@ -42,25 +42,25 @@ const GenericPropContainerImpl* AtomBase::GpImpl() const const String& AtomBase::GetName() const { this->CheckValidity(); - return impl_->GetName(); + return impl_->Name(); } void AtomBase::SetName(const String& atom_name) { this->CheckValidity(); - return impl_->SetName(atom_name); + impl_->Name()=atom_name; } const geom::Vec3& AtomBase::GetPos() const { this->CheckValidity(); - return impl_->GetPos(); + return impl_->TransformedPos(); } const geom::Vec3& AtomBase::GetOriginalPos() const { this->CheckValidity(); - return impl_->GetOriginalPos(); + return impl_->OriginalPos(); } geom::Vec3 AtomBase::GetAltPos(const String& alt_group) const diff --git a/modules/mol/base/src/editor_base.cc b/modules/mol/base/src/editor_base.cc index 01eaf549d..5ab540f31 100644 --- a/modules/mol/base/src/editor_base.cc +++ b/modules/mol/base/src/editor_base.cc @@ -143,7 +143,7 @@ void EditorBase::ReorderAllResidues() void EditorBase::RenameAtom(AtomHandle atom, const String& new_name) { CheckHandleValidity(atom); - atom.Impl()->SetName(new_name); + atom.Impl()->Name()=new_name; } BondHandle EditorBase::Connect(const AtomHandle& first, diff --git a/modules/mol/base/src/impl/atom_impl.cc b/modules/mol/base/src/impl/atom_impl.cc index a04f4ce1e..73123d538 100644 --- a/modules/mol/base/src/impl/atom_impl.cc +++ b/modules/mol/base/src/impl/atom_impl.cc @@ -95,11 +95,11 @@ void AtomImpl::TraceDirectionality(FragmentImplP frag, ConnectorImplP conn, #if !defined(NDEBUG) if (conn->GetFirst()==shared_from_this()) { LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber() - << "." << GetName() << " [" << conn->GetSecond()->GetQualifiedName() + << "." << Name() << " [" << conn->GetSecond()->GetQualifiedName() << " ]"); } else { LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber() - << "." << GetName() << " [" << conn->GetFirst()->GetQualifiedName() + << "." << Name() << " [" << conn->GetFirst()->GetQualifiedName() << " ]"); } @@ -107,7 +107,7 @@ void AtomImpl::TraceDirectionality(FragmentImplP frag, ConnectorImplP conn, #endif } else { LOG_TRACE("dir:" << String(n,' ') << " atom " << res_.lock()->GetNumber() - << "." << GetName() << " [ ]"); + << "." << Name() << " [ ]"); } // presence of a primary connector indicates ring closure @@ -196,7 +196,7 @@ void AtomImpl::UpdateFromXCS() // stack before calling UpdateFromICS() on the next atom. { // Derive direction and length of connector from atom positions. - geom::Vec3 global_d=((*i)->GetSecond()->GetOriginalPos()-this->GetOriginalPos()); + geom::Vec3 global_d=((*i)->GetSecond()->OriginalPos()-this->OriginalPos()); // Set direction and length of connector. Direction is relative to // local coordinate system of this atom. // Note the order of arguments for the matrix multiplication. This is the @@ -225,7 +225,7 @@ std::ostream& operator<<(std::ostream& o, const AtomImplPtr ap) { o << ap->GetResidue()->GetChain()->GetName() << "."; o << ap->GetResidue()->GetKey() << ap->GetResidue()->GetNumber() << "."; - o << ap->GetName(); + o << ap->Name(); return o; } @@ -258,7 +258,7 @@ ConnectorImplP GetConnector(const AtomImplPtr& a, const AtomImplPtr& b) { } String AtomImpl::GetQualifiedName() const { - return this->GetResidue()->GetQualifiedName()+"."+this->GetName(); + return this->GetResidue()->GetQualifiedName()+"."+this->Name(); } void AtomImpl::DeleteAllConnectors() { diff --git a/modules/mol/base/src/impl/atom_impl.hh b/modules/mol/base/src/impl/atom_impl.hh index 963f5ab60..b342c526a 100644 --- a/modules/mol/base/src/impl/atom_impl.hh +++ b/modules/mol/base/src/impl/atom_impl.hh @@ -56,21 +56,21 @@ public: ~AtomImpl(); void Apply(EntityVisitor& h); - const String& GetName() const {return name_;} - - void SetName(const String& atom_name) { - name_=atom_name; - } + // for efficiency reasons, the simple setter/getter methods are + // replaced by direct access - this is the impl layer after all + const String& Name() const {return name_;} + String& Name() {return name_;} - const geom::Vec3& GetPos() const {return tf_pos_;} + // DEPRECATED + const String& GetName() const {return name_;} - const geom::Vec3& GetOriginalPos() const {return pos_;} + const geom::Vec3& TransformedPos() const {return tf_pos_;} + geom::Vec3& TransformedPos() {return tf_pos_;} - void SetTransformedPos(const geom::Vec3& pos) { tf_pos_=pos; } + const geom::Vec3& OriginalPos() const {return pos_;} + geom::Vec3& OriginalPos() {return pos_;} - void SetOriginalPos(const geom::Vec3& pos) { pos_=pos; } - ResidueImplPtr GetResidue() const; void SetPrimaryConnector(const ConnectorImplP& bp) { diff --git a/modules/mol/base/src/impl/chain_impl.cc b/modules/mol/base/src/impl/chain_impl.cc index 456255c3c..ea0f4cce8 100644 --- a/modules/mol/base/src/impl/chain_impl.cc +++ b/modules/mol/base/src/impl/chain_impl.cc @@ -387,8 +387,8 @@ geom::AlignedCuboid ChainImpl::GetBounds() const ResidueImplPtr r=*i; for (AtomImplList::iterator j=r->GetAtomList().begin(); j!=r->GetAtomList().end(); ++j) { - mmin=geom::Min(mmin, (*j)->GetPos()); - mmax=geom::Max(mmax, (*j)->GetPos()); + mmin=geom::Min(mmin, (*j)->TransformedPos()); + mmax=geom::Max(mmax, (*j)->TransformedPos()); atoms=true; } } @@ -408,7 +408,7 @@ geom::Vec3 ChainImpl::GetCenterOfAtoms() const ResidueImplPtr r=*i; for (AtomImplList::iterator j=r->GetAtomList().begin(); j!=r->GetAtomList().end(); ++j) { - sum+=(*j)->GetPos(); + sum+=(*j)->TransformedPos(); } } sum/=this->GetAtomCount(); @@ -426,7 +426,7 @@ geom::Vec3 ChainImpl::GetCenterOfMass() const ResidueImplPtr r=*i; for (AtomImplList::iterator j=r->GetAtomList().begin(); j!=r->GetAtomList().end(); ++j) { - center+=(*j)->GetPos() * (*j)->GetMass(); + center+=(*j)->TransformedPos() * (*j)->GetMass(); } } center/=mass; diff --git a/modules/mol/base/src/impl/connector_impl.cc b/modules/mol/base/src/impl/connector_impl.cc index b244bc000..7598fdf6f 100644 --- a/modules/mol/base/src/impl/connector_impl.cc +++ b/modules/mol/base/src/impl/connector_impl.cc @@ -77,14 +77,14 @@ geom::Mat3 find_rotation(const geom::Vec3& d) { geom::Vec3 ConnectorImpl::GetPos() const { - return (GetFirst()->GetPos()+GetSecond()->GetPos())*0.5; + return (GetFirst()->TransformedPos()+GetSecond()->TransformedPos())*0.5; } Real ConnectorImpl::GetLength() const { Real length; if (this->GetFirst()->GetEntity()->HasICS()==false) { - length=geom::Distance(this->GetFirst()->GetOriginalPos(),this->GetSecond()->GetOriginalPos()); + length=geom::Distance(this->GetFirst()->OriginalPos(),this->GetSecond()->OriginalPos()); } else { length=len_; } @@ -110,8 +110,8 @@ bool ConnectorImpl::IsConnectorOf(const AtomImplPtr& a, geom::Vec3 ConnectorImpl::GetOriginalPos() const { - return (this->GetFirst()->GetOriginalPos()+ - this->GetSecond()->GetOriginalPos())*0.5; + return (this->GetFirst()->OriginalPos()+ + this->GetSecond()->OriginalPos())*0.5; } void ConnectorImpl::Switch() diff --git a/modules/mol/base/src/impl/dihedral.cc b/modules/mol/base/src/impl/dihedral.cc index 3580829b1..cac771039 100644 --- a/modules/mol/base/src/impl/dihedral.cc +++ b/modules/mol/base/src/impl/dihedral.cc @@ -46,8 +46,8 @@ Dihedral::Dihedral(const AtomImplList& atoms) Real Dihedral::GetAngleXCS() const { - return geom::DihedralAngle(atoms_[0]->GetPos(), atoms_[1]->GetPos(), - atoms_[2]->GetPos(), atoms_[3]->GetPos()); + return geom::DihedralAngle(atoms_[0]->TransformedPos(), atoms_[1]->TransformedPos(), + atoms_[2]->TransformedPos(), atoms_[3]->TransformedPos()); } Real Dihedral::GetAngleICS() const { diff --git a/modules/mol/base/src/impl/entity_impl.cc b/modules/mol/base/src/impl/entity_impl.cc index eb85cef79..622e06b1a 100644 --- a/modules/mol/base/src/impl/entity_impl.cc +++ b/modules/mol/base/src/impl/entity_impl.cc @@ -268,10 +268,10 @@ geom::AlignedCuboid EntityImpl::GetBounds() const if (this->GetAtomCount()>0) { geom::Vec3 mmin, mmax; AtomImplMap::const_iterator it=atom_map_.begin(); - mmin=mmax=it->second->GetPos(); + mmin=mmax=it->second->TransformedPos(); for (++it; it!=atom_map_.end();++it) { - mmin=geom::Min(mmin,it->second->GetPos()); - mmax=geom::Max(mmax,it->second->GetPos()); + mmin=geom::Min(mmin,it->second->TransformedPos()); + mmax=geom::Max(mmax,it->second->TransformedPos()); } return geom::AlignedCuboid(mmin, mmax); } else { @@ -285,7 +285,7 @@ geom::Vec3 EntityImpl::GetCenterOfAtoms() const { if (this->GetAtomCount()>0) { for (AtomImplMap::const_iterator it = atom_map_.begin(); it!=atom_map_.end();++it) { - center+=it->second->GetPos(); + center+=it->second->TransformedPos(); } center/=static_cast<Real>(atom_map_.size()); } @@ -297,7 +297,7 @@ geom::Vec3 EntityImpl::GetCenterOfMass() const { Real mass = this->GetMass(); if (this->GetAtomCount()>0 && mass>0) { for(AtomImplMap::const_iterator it = atom_map_.begin();it!=atom_map_.end();++it) { - center+=it->second->GetPos()*it->second->GetMass(); + center+=it->second->TransformedPos()*it->second->GetMass(); } center/=mass; } @@ -324,12 +324,12 @@ AtomImplPtr EntityImpl::CreateAtom(const ResidueImplPtr& rp, #else AtomImplPtr ap(new AtomImpl(shared_from_this(), rp, name, pos, ele, next_index_++)); #endif - if (identity_transf_ == false) { + if (!identity_transf_) { geom::Vec3 transformed_pos = geom::Vec3(transformation_matrix_*geom::Vec4(pos)); - ap->SetTransformedPos(transformed_pos); + ap->TransformedPos()=transformed_pos; atom_organizer_.Add(ap,transformed_pos); } else { - ap->SetTransformedPos(pos); + ap->TransformedPos()=pos; atom_organizer_.Add(ap,pos); } atom_map_.insert(AtomImplMap::value_type(ap.get(),ap)); @@ -496,8 +496,8 @@ Real EntityImpl::GetAngleXCS(const AtomImplPtr& a1, const AtomImplPtr& a2, ConnectorImplP c23=GetConnector(a2, a3); if (c12 && c12->GetFirst() && c12->GetSecond()) { if (c23 && c23->GetFirst() && c23->GetSecond()) { - return Angle(a2->GetPos()-a1->GetPos(), - a2->GetPos()-a3->GetPos()); + return Angle(a2->TransformedPos()-a1->TransformedPos(), + a2->TransformedPos()-a3->TransformedPos()); } else { AtomHandle ah2(a2), ah3(a3); throw NotConnectedError(ah2, ah3); @@ -763,12 +763,7 @@ void EntityImpl::SetTransform(const geom::Mat4 transfmat) { transformation_matrix_=transfmat; inverse_transformation_matrix_=Invert(transformation_matrix_); - geom::Mat4 identity = geom::Mat4(); - if (transformation_matrix_ == identity) { - identity_transf_ = true; - } else { - identity_transf_ = false; - } + identity_transf_ = (transformation_matrix_==geom::Mat4()); this->UpdateTransformedPos(); this->MarkXCSDirty(); } @@ -1047,7 +1042,7 @@ void EntityImpl::UpdateOrganizer() atom_organizer_.Clear(); for (AtomImplMap::const_iterator i=atom_map_.begin(), e=atom_map_.end(); i!=e; ++i) { - atom_organizer_.Add(i->second, i->second->GetPos()); + atom_organizer_.Add(i->second, i->second->TransformedPos()); } } @@ -1186,7 +1181,7 @@ void EntityImpl::RenameChain(ChainImplPtr chain, const String& new_name) void EntityImpl::UpdateTransformedPos(){ for(AtomImplMap::iterator it = atom_map_.begin();it!=atom_map_.end();++it) { - it->second->SetTransformedPos(geom::Vec3(transformation_matrix_*geom::Vec4(it->second->GetOriginalPos()))); + it->second->TransformedPos()=geom::Vec3(transformation_matrix_*geom::Vec4(it->second->OriginalPos())); } } diff --git a/modules/mol/base/src/impl/residue_impl.cc b/modules/mol/base/src/impl/residue_impl.cc index 3b8bd2ace..af7012825 100644 --- a/modules/mol/base/src/impl/residue_impl.cc +++ b/modules/mol/base/src/impl/residue_impl.cc @@ -60,8 +60,8 @@ AtomImplPtr ResidueImpl::InsertAtom(const String& name, AtomImplPtr ResidueImpl::InsertAtom(const AtomImplPtr& atom) { - AtomImplPtr dst_atom=this->InsertAtom(atom->GetName(), - atom->GetPos(), + AtomImplPtr dst_atom=this->InsertAtom(atom->Name(), + atom->TransformedPos(), atom->GetElement()); dst_atom->Assign(*atom.get()); @@ -166,12 +166,12 @@ AtomImplPtr ResidueImpl::GetCentralAtom() const if (chem_class_.IsNucleotideLinking()) { for (AtomImplList::const_iterator it=atom_list_.begin(); it!=atom_list_.end();++it) { - if((*it)->GetName()=="P") return *it; + if((*it)->Name()=="P") return *it; } } else if (chem_class_.IsPeptideLinking()) { for (AtomImplList::const_iterator it=atom_list_.begin(); it!=atom_list_.end();++it) { - if((*it)->GetName()=="CA") return *it; + if((*it)->Name()=="CA") return *it; } } @@ -217,14 +217,14 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const AtomImplPtr a1 = FindAtom("C"); AtomImplPtr a2 = FindAtom("O"); if(a1 && a2) { - nrvo = geom::Normalize(a2->GetPos()-a1->GetPos()); + nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos()); } else { a1 = FindAtom("CB"); a2 = FindAtom("CA"); if(a1 && a2) { - nrvo = geom::Normalize(a2->GetPos()-a1->GetPos()); + nrvo = geom::Normalize(a2->TransformedPos()-a1->TransformedPos()); } else { - geom::Vec3 v0=GetCentralAtom()->GetPos(); + geom::Vec3 v0=GetCentralAtom()->TransformedPos(); nrvo=geom::Cross(geom::Normalize(v0), geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1]))); LOG_VERBOSE("warning: could not find atoms for proper central normal calculation"); @@ -235,9 +235,9 @@ geom::Vec3 ResidueImpl::GetCentralNormal() const AtomImplPtr a2 = FindAtom("OP1"); AtomImplPtr a3 = FindAtom("OP2"); if(a1 && a2 && a3) { - nrvo = geom::Normalize(a1->GetPos()-(a2->GetPos()+a3->GetPos())*.5); + nrvo = geom::Normalize(a1->TransformedPos()-(a2->TransformedPos()+a3->TransformedPos())*.5); } else { - geom::Vec3 v0=GetCentralAtom()->GetPos(); + geom::Vec3 v0=GetCentralAtom()->TransformedPos(); nrvo=geom::Cross(geom::Normalize(v0), geom::Normalize(geom::Vec3(-v0[2],v0[0],v0[1]))); LOG_VERBOSE("warning: could not find atoms for proper central normal calculation"); @@ -251,7 +251,7 @@ AtomImplPtr ResidueImpl::FindAtom(const String& aname) const { for (AtomImplList::const_iterator it=atom_list_.begin(); it!=atom_list_.end();++it) { - if ((*it)->GetName()==aname) { + if ((*it)->Name()==aname) { return *it; } } @@ -366,18 +366,25 @@ void ResidueImpl::DeleteAtom(const AtomImplPtr& atom) { } } +namespace { + struct aname_matcher { + aname_matcher(const String& n): aname(n) {} + bool operator()(AtomImplPtr& a) {return a->Name()==aname;} + String aname; + }; +} + void ResidueImpl::DeleteAtoms(const String& atom_name) { AtomImplList::iterator i=atom_list_.begin(); for (; i!=atom_list_.end(); ++i) { - if ((*i)->GetName()==atom_name) { + if ((*i)->Name()==atom_name) { (*i)->DeleteAllTorsions(); (*i)->DeleteAllConnectors(); this->GetEntity()->DeleteAtom(*i); } } AtomImplList::iterator new_end; - new_end=std::remove_if(atom_list_.begin(), atom_list_.end(), - bind(&AtomImpl::GetName, _1)==atom_name); + new_end=std::remove_if(atom_list_.begin(), atom_list_.end(), aname_matcher(atom_name)); atom_list_.erase(new_end, atom_list_.end()); } @@ -452,10 +459,10 @@ geom::AlignedCuboid ResidueImpl::GetBounds() const if (atom_list_.size()>0) { AtomImplList::const_iterator i=atom_list_.begin(); - mmin=mmax=(*i)->GetPos(); + mmin=mmax=(*i)->TransformedPos(); for (++i; i!=atom_list_.end(); ++i) { - mmax=geom::Max(mmax,(*i)->GetPos()); - mmin=geom::Min(mmin,(*i)->GetPos()); + mmax=geom::Max(mmax,(*i)->TransformedPos()); + mmin=geom::Min(mmin,(*i)->TransformedPos()); } return geom::AlignedCuboid(mmin, mmax); } else { @@ -469,7 +476,7 @@ geom::Vec3 ResidueImpl::GetCenterOfAtoms() const if (!atom_list_.empty()) { for (AtomImplList::const_iterator i=atom_list_.begin(); i!=atom_list_.end(); ++i) { - sum+=(*i)->GetPos(); + sum+=(*i)->TransformedPos(); } sum/=atom_list_.size(); } @@ -483,7 +490,7 @@ geom::Vec3 ResidueImpl::GetCenterOfMass() const if (this->GetAtomCount() > 0 && mass > 0) { for (AtomImplList::const_iterator i=atom_list_.begin(); i!=atom_list_.end(); ++i) { - center+=(*i)->GetPos()*(*i)->GetMass(); + center+=(*i)->TransformedPos()*(*i)->GetMass(); } } return center/mass; @@ -565,7 +572,7 @@ bool ResidueImpl::SwitchAtomPos(const String& group) { for (; k!=gr.atoms.end(); ++k) { AtomGroupEntry& entry=*k; assert(!entry.atom.expired()); - entry.pos=entry.atom.lock()->GetOriginalPos(); + entry.pos=entry.atom.lock()->OriginalPos(); } } AtomGroup& agr=i->second; @@ -574,11 +581,11 @@ bool ResidueImpl::SwitchAtomPos(const String& group) { AtomGroupEntry& entry=*j; assert(!entry.atom.expired()); - entry.atom.lock()->SetOriginalPos(entry.pos); + entry.atom.lock()->OriginalPos()=entry.pos; EntityHandle ent = entry.atom.lock()->GetEntity(); geom::Mat4 transf_matrix = ent.GetTransformationMatrix(); geom::Vec3 transf_pos = geom::Vec3(transf_matrix*geom::Vec4(entry.pos)); - entry.atom.lock()->SetTransformedPos(transf_pos); + entry.atom.lock()->TransformedPos()=transf_pos; } curr_group_=group; return true; diff --git a/modules/mol/base/src/query_state.cc b/modules/mol/base/src/query_state.cc index 2f04aeecc..4e5c86403 100644 --- a/modules/mol/base/src/query_state.cc +++ b/modules/mol/base/src/query_state.cc @@ -318,22 +318,22 @@ boost::logic::tribool QueryState::EvalAtom(const impl::AtomImplPtr& a) { int int_value; switch (ss.sel_id) { case Prop::ANAME: - str_value = a->GetName(); + str_value = a->Name(); s_[*i] = cmp_string(ss.comp_op,str_value, boost::get<String>(ss.param)); break; case Prop::AX: - float_value=(a->GetPos())[0]; + float_value=(a->TransformedPos())[0]; s_[*i]=cmp_num<Real>(ss.comp_op, float_value, boost::get<float>(ss.param)); break; case Prop::AY: - float_value=(a->GetPos())[1]; + float_value=(a->TransformedPos())[1]; s_[*i]=cmp_num<Real>(ss.comp_op, float_value, boost::get<float>(ss.param)); break; case Prop::AZ: - float_value=(a->GetPos())[2]; + float_value=(a->TransformedPos())[2]; s_[*i]=cmp_num<Real>(ss.comp_op, float_value, boost::get<float>(ss.param)); break; @@ -353,7 +353,7 @@ boost::logic::tribool QueryState::EvalAtom(const impl::AtomImplPtr& a) { boost::get<float>(ss.param)); break; case Prop::WITHIN: - s_[*i]= this->do_within(a->GetPos(), + s_[*i]= this->do_within(a->TransformedPos(), boost::get<WithinParam>(ss.param), ss.comp_op); break; diff --git a/modules/mol/base/src/xcs_editor.cc b/modules/mol/base/src/xcs_editor.cc index 88d8bbcd4..41d6be611 100644 --- a/modules/mol/base/src/xcs_editor.cc +++ b/modules/mol/base/src/xcs_editor.cc @@ -67,33 +67,117 @@ XCSEditor& XCSEditor::operator=(const XCSEditor& rhs) return *this; } -void XCSEditor::SetAtomPos(const AtomHandle& atom, - const geom::Vec3& position) +void XCSEditor::SetAtomTransformedPos(const AtomHandle& atom, + const geom::Vec3& position) { CheckHandleValidity(atom); - atom.Impl()->SetTransformedPos(position); - geom::Mat4 inv_transformation_matrix = ent_.Impl()->GetInvTransfMatrix(); - geom::Vec3 original_pos = geom::Vec3(inv_transformation_matrix*geom::Vec4(position)); - atom.Impl()->SetOriginalPos(original_pos); + atom.Impl()->TransformedPos()=position; + if(ent_.Impl()->IsTransfIdentity()) { + atom.Impl()->OriginalPos()=position; + } else { + atom.Impl()->OriginalPos() = geom::Vec3(ent_.Impl()->GetInvTransfMatrix()*geom::Vec4(position)); + } ent_.Impl()->MarkICSDirty(); ent_.Impl()->MarkOrganizerDirty(); this->Update(); } +namespace { + template<typename T> + void set_transformed_pos(impl::EntityImpl* ent, const AtomHandleList& alist, T *positions) + { + bool has_tf=ent->IsTransfIdentity(); + for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait) { + ait->Impl()->TransformedPos()[0]=static_cast<Real>(positions[0]); + ait->Impl()->TransformedPos()[1]=static_cast<Real>(positions[1]); + ait->Impl()->TransformedPos()[2]=static_cast<Real>(positions[2]); + positions+=3; + if(has_tf) { + ait->Impl()->OriginalPos()=ait->Impl()->TransformedPos(); + } else { + ait->Impl()->OriginalPos() = geom::Vec3(ent->GetInvTransfMatrix()*geom::Vec4(ait->Impl()->TransformedPos())); + } + } + ent->MarkICSDirty(); + ent->MarkOrganizerDirty(); + } +} // anon ns + +void XCSEditor::SetAtomTransformedPos(const AtomHandleList& alist, float *positions) +{ + set_transformed_pos<float>(ent_.Impl().get(),alist,positions); + this->Update(); +} + +void XCSEditor::SetAtomTransformedPos(const AtomHandleList& alist, double *positions) +{ + set_transformed_pos<double>(ent_.Impl().get(),alist,positions); + this->Update(); +} void XCSEditor::SetAtomOriginalPos(const AtomHandle& atom, - const geom::Vec3& position) + const geom::Vec3& position) { CheckHandleValidity(atom); - atom.Impl()->SetOriginalPos(position); - geom::Mat4 transformation_matrix = atom.GetEntity().GetTransformationMatrix(); - geom::Vec3 transformed_pos = geom::Vec3(transformation_matrix*geom::Vec4(position)); - atom.Impl()->SetTransformedPos(transformed_pos); + atom.Impl()->OriginalPos()=position; + if(ent_.Impl()->IsTransfIdentity()) { + atom.Impl()->TransformedPos()=position; + } else { + atom.Impl()->TransformedPos() = geom::Vec3(ent_.Impl()->GetTransfMatrix()*geom::Vec4(position)); + } ent_.Impl()->MarkICSDirty(); ent_.Impl()->MarkOrganizerDirty(); this->Update(); } +namespace { + template<typename T> + void set_original_pos(impl::EntityImpl* ent, const AtomHandleList& alist, T *positions) + { + bool has_tf=ent->IsTransfIdentity(); + for(AtomHandleList::const_iterator ait=alist.begin();ait!=alist.end();++ait) { + ait->Impl()->OriginalPos()[0]=static_cast<Real>(positions[0]); + ait->Impl()->OriginalPos()[1]=static_cast<Real>(positions[1]); + ait->Impl()->OriginalPos()[2]=static_cast<Real>(positions[2]); + positions+=3; + if(has_tf) { + ait->Impl()->TransformedPos()=ait->Impl()->OriginalPos(); + } else { + ait->Impl()->TransformedPos() = geom::Vec3(ent->GetTransfMatrix()*geom::Vec4(ait->Impl()->OriginalPos())); + } + } + ent->MarkICSDirty(); + ent->MarkOrganizerDirty(); + } +} // anon ns + +void XCSEditor::SetAtomOriginalPos(const AtomHandleList& alist, float *positions) +{ + set_original_pos<float>(ent_.Impl().get(),alist,positions); + this->Update(); +} + +void XCSEditor::SetAtomOriginalPos(const AtomHandleList& alist, double *positions) +{ + set_original_pos<double>(ent_.Impl().get(),alist,positions); + this->Update(); +} + +void XCSEditor::SetAtomPos(const AtomHandle& atom, const geom::Vec3& position) +{ + this->SetAtomTransformedPos(atom,position); +} + +void XCSEditor::SetAtomPos(const AtomHandleList& alist, float *positions) +{ + this->SetAtomTransformedPos(alist,positions); +} + +void XCSEditor::SetAtomPos(const AtomHandleList& alist, double *positions) +{ + this->SetAtomTransformedPos(alist,positions); +} + void XCSEditor::ApplyTransform(const geom::Mat4& transform) { ent_.Impl()->ApplyTransform(transform); diff --git a/modules/mol/base/src/xcs_editor.hh b/modules/mol/base/src/xcs_editor.hh index 4b89083c6..73d96a652 100644 --- a/modules/mol/base/src/xcs_editor.hh +++ b/modules/mol/base/src/xcs_editor.hh @@ -51,11 +51,47 @@ public: void SetAtomOriginalPos(const AtomHandle& atom, const geom::Vec3& position); + /// \brief numpy float interface + /// + /// the passed in float array must have a length of 3*alist.size() + void SetAtomOriginalPos(const AtomHandleList& alist, + float *positions); + + /// \brief numpy double interface + /// + /// the passed in double array must have a length of 3*alist.size() + void SetAtomOriginalPos(const AtomHandleList& alist, + double *positions); + /// \brief set transformed position of atom /// /// This function also updates the original coordinates + void SetAtomTransformedPos(const AtomHandle& atom, + const geom::Vec3& position); + + /// \brief numpy float interface + /// + /// the passed in float array must have a length of 3*alist.size() + void SetAtomTransformedPos(const AtomHandleList& alist, + float *positions); + + /// \brief numpy double interface + /// + /// the passed in double array must have a length of 3*alist.size() + void SetAtomTransformedPos(const AtomHandleList& alist, + double *positions); + + /// \brief same as SetAtomTransformedPos(AtomHandle, geom::Vec3) void SetAtomPos(const AtomHandle& atom, - const geom::Vec3& position); + const geom::Vec3& position); + + /// \brief same as SetTransformedPos(AtomHandleList,float*) + void SetAtomPos(const AtomHandleList& alist, + float *positions); + + /// \brief same as SetTransformedPos(AtomHandleList,double*) + void SetAtomPos(const AtomHandleList& alist, + double *positions); /// \brief apply additional transformation to all atoms /// @@ -69,6 +105,7 @@ public: /// \brief immediately update internal coordinate system void UpdateICS(); + protected: XCSEditor(const EntityHandle& ent, EditMode mode); diff --git a/modules/mol/base/tests/test_numpy.py b/modules/mol/base/tests/test_numpy.py index 8dc14835a..cad6f764f 100644 --- a/modules/mol/base/tests/test_numpy.py +++ b/modules/mol/base/tests/test_numpy.py @@ -4,6 +4,9 @@ sys.path.append("../../../../stage/lib64/openstructure") from ost import * import numpy +def v2v(v): + return geom.Vec3(float(v[0]),float(v[1]),float(v[2])) + class TestNumpy(unittest.TestCase): def setUp(self): pass @@ -14,16 +17,20 @@ class TestNumpy(unittest.TestCase): ch=ed.InsertChain("X") re=ed.AppendResidue(ch,"ALA") a0=ed.InsertAtom(re,"A",geom.Vec3(0,0,0)) + self.assertTrue(a0.GetIndex()==0) a1=ed.InsertAtom(re,"B",geom.Vec3(1,0,0)) + self.assertTrue(a1.GetIndex()==1) a2=ed.InsertAtom(re,"C",geom.Vec3(2,0,0)) + self.assertTrue(a2.GetIndex()==2) a3=ed.InsertAtom(re,"D",geom.Vec3(3,0,0)) + self.assertTrue(a3.GetIndex()==3) self.assertTrue(geom.Distance(a0.pos,geom.Vec3(0,0,0))<1e-10) self.assertTrue(geom.Distance(a1.pos,geom.Vec3(1,0,0))<1e-10) self.assertTrue(geom.Distance(a2.pos,geom.Vec3(2,0,0))<1e-10) self.assertTrue(geom.Distance(a3.pos,geom.Vec3(3,0,0))<1e-10) - entity.positions=numpy.array([[0,1,0],[0,2,0],[0,3,0],[0,4,0]]) + ed.SetAtomTransformedPos(entity.GetAtomList(),numpy.array([[0,1,0],[0,2,0],[0,3,0],[0,4,0]], dtype=numpy.float32)) self.assertTrue(geom.Distance(a0.pos,geom.Vec3(0,1,0))<1e-10) self.assertTrue(geom.Distance(a1.pos,geom.Vec3(0,2,0))<1e-10) @@ -32,10 +39,10 @@ class TestNumpy(unittest.TestCase): na=entity.positions - self.assertTrue(geom.Distance(geom.Vec3(na[0][0],na[0][1],na[0][2]),geom.Vec3(0,1,0))<1e-10) - self.assertTrue(geom.Distance(geom.Vec3(na[1][0],na[1][1],na[1][2]),geom.Vec3(0,2,0))<1e-10) - self.assertTrue(geom.Distance(geom.Vec3(na[2][0],na[2][1],na[2][2]),geom.Vec3(0,3,0))<1e-10) - self.assertTrue(geom.Distance(geom.Vec3(na[3][0],na[3][1],na[3][2]),geom.Vec3(0,4,0))<1e-10) + self.assertTrue(geom.Distance(v2v(na[0]),geom.Vec3(0,1,0))<1e-10) + self.assertTrue(geom.Distance(v2v(na[1]),geom.Vec3(0,2,0))<1e-10) + self.assertTrue(geom.Distance(v2v(na[2]),geom.Vec3(0,3,0))<1e-10) + self.assertTrue(geom.Distance(v2v(na[3]),geom.Vec3(0,4,0))<1e-10) if __name__== '__main__': unittest.main() -- GitLab