backup -- met nieuwe Cif++ optimalisaties

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@341 a1961a4f-ab94-4bcc-80e8-33b5a54de466

include/cif++/Cif++.h

@@ -611,6 +611,7 @@ struct ConditionImpl
 {
 	virtual ~ConditionImpl() {}
 	
+	virtual void prepare(const Category& c) {}
 	virtual bool test(const Category& c, const Row& r) const = 0;
 	virtual std::string str() const = 0;
 };
@@ -638,9 +639,16 @@ struct Condition
 		delete mImpl;
 	}
 	
+	void prepare(const Category& c)
+	{
+		mImpl->prepare(c);
+		mPrepared = true;
+	}
+	
 	bool operator()(const Category& c, const Row& r) const
 	{
 		assert(mImpl);
+		assert(mPrepared);
 		return mImpl->test(c, r);
 	}
 	
@@ -650,6 +658,7 @@ struct Condition
 	}
 
 	detail::ConditionImpl*	mImpl;
+	bool mPrepared = false;
 };
 
 namespace detail
@@ -660,9 +669,11 @@ struct KeyIsEmptyConditionImpl : public ConditionImpl
 	KeyIsEmptyConditionImpl(const string& ItemTag)
 		: mItemTag(ItemTag) {}
 
+	virtual void prepare(const Category& c);
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
-		return r[mItemTag].empty();
+		return r[mItemIx].empty();
 	}
 	
 	virtual std::string str() const
@@ -671,6 +682,7 @@ struct KeyIsEmptyConditionImpl : public ConditionImpl
 	}
 	
 	string mItemTag;
+	size_t mItemIx;
 };
 
 template<typename T>
@@ -681,9 +693,11 @@ struct KeyIsConditionImpl : public ConditionImpl
 	KeyIsConditionImpl(const string& ItemTag, const valueType& value)
 		: mItemTag(ItemTag), mValue(value) {}
 	
+	virtual void prepare(const Category& c);
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
-		return r[mItemTag].template compare<valueType>(mValue) == 0;
+		return r[mItemIx].template compare<valueType>(mValue) == 0;
 	}
 	
 	virtual std::string str() const
@@ -692,6 +706,7 @@ struct KeyIsConditionImpl : public ConditionImpl
 	}
 	
 	string mItemTag;
+	size_t mItemIx;
 	valueType mValue;
 };
 
@@ -703,9 +718,11 @@ struct KeyIsNotConditionImpl : public ConditionImpl
 	KeyIsNotConditionImpl(const string& ItemTag, const valueType& value)
 		: mItemTag(ItemTag), mValue(value) {}
 	
+	virtual void prepare(const Category& c);
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
-		return r[mItemTag].template compare<valueType>(mValue) != 0;
+		return r[mItemIx].template compare<valueType>(mValue) != 0;
 	}
 	
 	virtual std::string str() const
@@ -714,6 +731,7 @@ struct KeyIsNotConditionImpl : public ConditionImpl
 	}
 	
 	string mItemTag;
+	size_t mItemIx;
 	valueType mValue;
 };
 
@@ -723,6 +741,8 @@ struct KeyCompareConditionImpl : public ConditionImpl
 	KeyCompareConditionImpl(const string& ItemTag, COMP&& comp)
 		: mItemTag(ItemTag), mComp(std::move(comp)) {}
 	
+	virtual void prepare(const Category& c);
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
 		return mComp(c, r);
@@ -734,6 +754,7 @@ struct KeyCompareConditionImpl : public ConditionImpl
 	}
 	
 	string mItemTag;
+	size_t mItemIx;
 	COMP mComp;
 };
 
@@ -742,9 +763,11 @@ struct KeyMatchesConditionImpl : public ConditionImpl
 	KeyMatchesConditionImpl(const string& ItemTag, const std::regex& rx)
 		: mItemTag(ItemTag), mRx(rx) {}
 	
+	virtual void prepare(const Category& c);
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
-		return std::regex_match(r[mItemTag].as<string>(), mRx);
+		return std::regex_match(r[mItemIx].as<string>(), mRx);
 	}
 	
 	virtual std::string str() const
@@ -753,6 +776,7 @@ struct KeyMatchesConditionImpl : public ConditionImpl
 	}
 	
 	string mItemTag;
+	size_t mItemIx;
 	std::regex mRx;
 };
 
@@ -804,6 +828,12 @@ struct andConditionImpl : public ConditionImpl
 		delete mB;
 	}
 	
+	virtual void prepare(const Category& c)
+	{
+		mA->prepare(c);
+		mB->prepare(c);
+	}
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
 		return mA->test(c, r) and mB->test(c, r);
@@ -833,6 +863,12 @@ struct orConditionImpl : public ConditionImpl
 		delete mB;
 	}
 	
+	virtual void prepare(const Category& c)
+	{
+		mA->prepare(c);
+		mB->prepare(c);
+	}
+	
 	virtual bool test(const Category& c, const Row& r) const
 	{
 		return mA->test(c, r) or mB->test(c, r);
@@ -1218,22 +1254,40 @@ inline bool anyMatchesConditionImpl::test(const Category& c, const Row& r) const
 	
 }
 
-}
-
-namespace std
-{
+// these should be here, as I learned today
 
-template<>
 inline void swap(cif::Row& a, cif::Row& b)
 {
 	a.swap(b);
 }
 
-template<>
 inline void swap(cif::detail::ItemReference& a, cif::detail::ItemReference& b)
 {
 	a.swap(b);
 }
 
+namespace detail
+{
+
+template<typename T>
+void KeyIsConditionImpl<T>::prepare(const Category& c)
+{
+	mItemIx = c.getColumnIndex(mItemTag);
+}
+
+template<typename T>
+void KeyIsNotConditionImpl<T>::prepare(const Category& c)
+{
+	mItemIx = c.getColumnIndex(mItemTag);
+}
+
+template<typename T>
+void KeyCompareConditionImpl<T>::prepare(const Category& c)
+{
+	mItemIx = c.getColumnIndex(mItemTag);
+}
+
+}
+
 }
 

include/cif++/Structure.h

@@ -152,10 +152,20 @@ class Atom
 			labelAtomId() == "C" or labelAtomId() == "CA";
 	}
 	
+	void swap(Atom& b)
+	{
+		std::swap(mImpl, b.mImpl);
+	}
+
   private:
  	struct AtomImpl*			mImpl;
 };
 
+inline void swap(mmcif::Atom& a, mmcif::Atom& b)
+{
+	a.swap(b);
+}
+
 inline double Distance(const Atom& a, const Atom& b)
 {
 	return Distance(a.location(), b.location());
@@ -187,7 +197,7 @@ class Residue
 		, mAsymID(asymID), mAltID(altID), mSeqID(seqID) {}
 	
 	const Compound&		compound() const;
-	AtomView			atoms() const;
+	const AtomView&		atoms() const;
 	
 	Atom				atomByID(const std::string& atomID) const;
 
@@ -385,7 +395,7 @@ class Structure
 
 	File& getFile() const;
 
-	AtomView atoms() const;
+	const AtomView& atoms() const;
 	AtomView waters() const;
 	
 	std::vector<Polymer> polymers() const;

src/Cif++.cpp

@@ -485,6 +485,24 @@ void Datablock::write(ostream& os, const vector<string>& order)
 
 // --------------------------------------------------------------------
 //
+
+namespace detail
+{
+	
+void KeyIsEmptyConditionImpl::prepare(const Category& c)
+{
+	mItemIx = c.getColumnIndex(mItemTag);
+}
+
+void KeyMatchesConditionImpl::prepare(const Category& c)
+{
+	mItemIx = c.getColumnIndex(mItemTag);
+}
+
+}
+
+// --------------------------------------------------------------------
+//
 //	class to compare two rows based on their keys.
 
 class RowComparator
@@ -1229,6 +1247,8 @@ Row Category::operator[](Condition&& cond)
 {
 	Row result;
 
+	cond.prepare(*this);
+	
 	for (auto r: *this)
 	{
 		if (cond(*this, r))
@@ -1244,6 +1264,9 @@ Row Category::operator[](Condition&& cond)
 RowSet Category::find(Condition&& cond)
 {
 	RowSet result(*this);
+	
+	cond.prepare(*this);
+	
 	for (auto r: *this)
 	{
 		if (cond(*this, r))
@@ -1256,6 +1279,8 @@ bool Category::exists(Condition&& cond)
 {
 	bool result = false;
 
+	cond.prepare(*this);
+	
 	for (auto r: *this)
 	{
 		if (cond(*this, r))
@@ -1380,6 +1405,8 @@ void Category::erase(Condition&& cond)
 {
 	RowSet remove(*this);
 	
+	cond.prepare(*this);
+
 	for (auto r: *this)
 	{
 		if (cond(*this, r))
@@ -2201,8 +2228,8 @@ File::File(boost::filesystem::path p, bool validate)
 File::File(File&& rhs)
 	: mHead(nullptr), mValidator(nullptr)
 {
-	swap(mHead, rhs.mHead);
-	swap(mValidator, rhs.mValidator);
+	std::swap(mHead, rhs.mHead);
+	std::swap(mValidator, rhs.mValidator);
 }
 
 File::~File()

src/Compound.cpp

@@ -1118,7 +1118,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(fs::path file, CompoundFactoryImpl* nex
 //{
 //	const char* clibdMon = getenv("CLIBD_MON");
 //	if (clibdMon == nullptr)
-//		throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
+//		throw runtime_error("Cannot locate peptide list, please source the CCP4 environment");
 //	mClibdMon = clibdMon;
 //}
 //
@@ -1305,7 +1305,7 @@ CompoundFactory::CompoundFactory()
 {
 	const char* clibdMon = getenv("CLIBD_MON");
 	if (clibdMon == nullptr)
-		throw runtime_error("Cannot locate peptide list, please souce the CCP4 environment");
+		throw runtime_error("Cannot locate peptide list, please source the CCP4 environment");
 	
 	fs::path db = fs::path(clibdMon) / "list" / "mon_lib_list.cif";
 

src/DistanceMap.cpp

@@ -253,7 +253,7 @@ float DistanceMap::operator()(const Atom& a, const Atom& b) const
 	}
 	
 	if (ixb < ixa)
-		swap(ixa, ixb);
+		std::swap(ixa, ixb);
 	
 	tuple<size_t,size_t> k{ ixa, ixb };
 

src/Secondary.cpp

@@ -8,6 +8,8 @@
 #include "cif++/Config.h"
 
 #include <numeric>
+#include <chrono>
+#include <iomanip>
 
 #include <boost/algorithm/string.hpp>
 
@@ -775,21 +777,63 @@ struct DSSPImpl
 	vector<Res>			mResidues;
 };
 
+ostream& operator<<(ostream& os, const chrono::duration<double>& t)
+{
+	uint64 s = static_cast<uint64>(trunc(t.count()));
+	if (s > 24 * 60 * 60)
+	{
+		uint32 days = s / (24 * 60 * 60);
+		os << days << "d ";
+		s %= 24 * 60 * 60;
+	}
+	
+	if (s > 60 * 60)
+	{
+		uint32 hours = s / (60 * 60);
+		os << hours << "h ";
+		s %= 60 * 60;
+	}
+	
+	if (s > 60)
+	{
+		uint32 minutes = s / 60;
+		os << minutes << "m ";
+		s %= 60;
+	}
+	
+	double ss = s + 1e-6 * (t.count() - s);
+	
+	os << fixed << setprecision(1) << ss << 's';
+
+	return os;
+}
+
+
 DSSPImpl::DSSPImpl(const Structure& s)
 	: mStructure(s)
 	, mPolymers(mStructure.polymers())
 {
+	if (VERBOSE)
+		cerr << "Calculating DSSP ";
+	
 	size_t nRes = accumulate(mPolymers.begin(), mPolymers.end(),
 		0.0, [](double s, auto& p) { return s + p.size(); });
 	
 	mResidues.reserve(nRes);
 
+cerr << "starting to copy atoms for DSSP" << endl;
+auto start = std::chrono::system_clock::now();
+	
 	for (auto& p: mPolymers)
 	{
 		for (auto m: p)
 			mResidues.emplace_back(move(m));
 	}
 
+auto end = std::chrono::system_clock::now();
+chrono::duration<double> diff = end - start;
+cerr << "Copying atoms took " << diff << " seconds" << endl;
+	
 	for (size_t i = 0; i + 1 < mResidues.size(); ++i)
 	{
 		mResidues[i].mNext = &mResidues[i + 1];
@@ -798,11 +842,18 @@ DSSPImpl::DSSPImpl(const Structure& s)
 		mResidues[i + 1].assignHydrogen();
 	}
 	
+	if (VERBOSE) cerr << ".";
 	CalculateHBondEnergies(mResidues);
+
+	if (VERBOSE) cerr << ".";
 	CalculateBetaSheets(mResidues);
+
+	if (VERBOSE) cerr << ".";
 	CalculateAlphaHelices(mResidues);
 	
-	if (VERBOSE)
+	if (VERBOSE) cerr << endl;
+
+	if (VERBOSE > 1)
 	{
 		for (auto& r: mResidues)
 		{
@@ -823,7 +874,7 @@ DSSPImpl::DSSPImpl(const Structure& s)
 			
 			string id = m.asymID() + ':' + to_string(m.seqID()) + '/' + m.compoundID();
 			
-			cout << id << string(12 - id.length(), ' ')
+			cerr << id << string(12 - id.length(), ' ')
 				 << char(r.mSecondaryStructure) << ' '
 				 << helix
 				 << endl;

src/Structure.cpp

@@ -631,13 +631,14 @@ const Compound& Residue::compound() const
 	return *result;
 }
 
-AtomView Residue::atoms() const
+const AtomView& Residue::atoms() const
 {
 	if (mStructure == nullptr)
 		throw runtime_error("Invalid Residue object");
 
 	if (mAtoms.empty())
 	{
+#if 0
 		for (auto& a: mStructure->atoms())
 		{
 			if (mSeqID > 0 and a.labelSeqId() != mSeqID)
@@ -652,28 +653,22 @@ AtomView Residue::atoms() const
 			
 			mAtoms.push_back(a);
 		}
-//
-//		auto& atomSites = mStructure->category("atom_site");
-//		
-//		auto query = cif::Key("label_asym_id") == mAsymID and cif::Key("label_comp_id") == mCompoundID;
-//		
-//		if (not mSeqID.empty())
-//			query = move(query) and cif::Key("label_seq_id") == mSeqID;
-//		
-//		if (not mAltID.empty())
-//			query = move(query) and (cif::Key("label_alt_id") == cif::Empty() or cif::Key("label_alt_id") == mAltID);
-//	
-//		auto cifAtoms = atomSites.find(move(query));
-//
-//		set<string> ids;
-//		for (auto cifAtom: cifAtoms)
-//			ids.insert(cifAtom["id"].as<string>());
-//		
-//		for (auto& a: mStructure->atoms())
-//		{
-//			if (ids.count(a.id()))
-//				mAtoms.push_back(a);
-//		}
+#else
+		auto& atomSites = mStructure->category("atom_site");
+		
+		auto query = cif::Key("label_asym_id") == mAsymID and cif::Key("label_comp_id") == mCompoundID;
+		
+		if (mSeqID != 0)
+			query = move(query) and cif::Key("label_seq_id") == mSeqID;
+		
+		if (not mAltID.empty())
+			query = move(query) and (cif::Key("label_alt_id") == cif::Empty() or cif::Key("label_alt_id") == mAltID);
+	
+		auto cifAtoms = atomSites.find(move(query));
+
+		for (auto cifAtom: cifAtoms)
+			mAtoms.push_back(mStructure->getAtomById(cifAtom["id"].as<string>()));
+#endif
 	}
 	
 	return mAtoms;
@@ -785,12 +780,22 @@ Monomer::Monomer(const Polymer& polymer, uint32 index, int seqID, const string&
 float Monomer::phi() const
 {
 	float result = 360;
+
+	try
+	{
 		if (mIndex > 0)
 		{
 			Monomer prev = mPolymer->operator[](mIndex - 1);
 			if (prev.mSeqID + 1 == mSeqID)
 				result = DihedralAngle(prev.C().location(), N().location(), CAlpha().location(), C().location()); 
 		}
+	}
+	catch (const exception& ex)
+	{
+		if (VERBOSE)
+			cerr << ex.what() << endl;
+	}
+
 
 	return result;
 }
@@ -798,12 +803,21 @@ float Monomer::phi() const
 float Monomer::psi() const
 {
 	float result = 360;
+	
+	try
+	{
 		if (mIndex + 1 < mPolymer->size())
 		{
 			Monomer next = mPolymer->operator[](mIndex + 1);
 			if (mSeqID + 1 == next.mSeqID)
 				result = DihedralAngle(N().location(), CAlpha().location(), C().location(), next.N().location()); 
 		}
+	}
+	catch (const exception& ex)
+	{
+		if (VERBOSE)
+			cerr << ex.what() << endl;
+	}
 
 	return result;
 }
@@ -812,6 +826,8 @@ float Monomer::alpha() const
 {
 	float result = 360;
 
+	try
+	{
 		if (mIndex > 1 and mIndex + 2 < mPolymer->size())
 		{
 			Monomer prev = mPolymer->operator[](mIndex - 1);
@@ -820,6 +836,12 @@ float Monomer::alpha() const
 			
 			result = DihedralAngle(prev.CAlpha().location(), CAlpha().location(), next.CAlpha().location(), nextNext.CAlpha().location());
 		}
+	}
+	catch (const exception& ex)
+	{
+		if (VERBOSE)
+			cerr << ex.what() << endl;
+	}
 
 	return result;
 }
@@ -1162,6 +1184,8 @@ struct StructureImpl
 			if (modelNr.empty() or modelNr.as<uint32>() == mModelNr)
 				mAtoms.emplace_back(new AtomImpl(f, a["id"].as<string>(), a));
 		}
+		
+		sort(mAtoms.begin(), mAtoms.end(), [](auto& a, auto& b) { return a.id() < b.id(); });
 	}
 	
 	StructureImpl(const StructureImpl& si)
@@ -1268,11 +1292,11 @@ void StructureImpl::swapAtoms(Atom& a1, Atom& a2)
 	
 	auto l1 = r1.front()["label_atom_id"];
 	auto l2 = r2.front()["label_atom_id"];
-	swap(l1, l2);
+	std::swap(l1, l2);
 	
 	auto l3 = r1.front()["auth_atom_id"];
 	auto l4 = r2.front()["auth_atom_id"];
-	swap(l3, l4);
+	std::swap(l3, l4);
 }
 
 void StructureImpl::moveAtom(Atom& a, Point p)
@@ -1346,7 +1370,7 @@ Structure::Structure(const Structure& rhs)
 {
 }
 
-AtomView Structure::atoms() const
+const AtomView& Structure::atoms() const
 {
 	return mImpl->mAtoms;
 }
@@ -1420,13 +1444,13 @@ vector<Residue> Structure::nonPolymers() const
 
 Atom Structure::getAtomById(string id) const
 {
-	for (auto& a: mImpl->mAtoms)
-	{
-		if (a.id() == id)
-			return a;
-	}
+	auto i = lower_bound(mImpl->mAtoms.begin(), mImpl->mAtoms.end(),
+		id, [](auto& a, auto& b) { return a.id() < b; });
 
+	if (i == mImpl->mAtoms.end() or i->id() != id)
 		throw out_of_range("Could not find atom with id " + id);
+
+	return *i;
 }
 
 File& Structure::getFile() const