Skip to content

L
libcifpp
Commit 065ba540 by maarten
Options

faster distmap 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@257 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent e19fc4f7
Showing with 99 additions and 33 deletions
include/cif++/BondMap.h
......@@ -21,7 +21,7 @@ class BondMap
private:
uint32 dim;
size_t dim;
std::vector<bool> bond;
std::unordered_map<std::string,size_t> index;
};
......
include/cif++/DistanceMap.h
......@@ -4,6 +4,8 @@
#include <unordered_map>
//#include <boost/hash/hash.hpp>
#include <clipper/clipper.h>
#include "cif++/Structure.h"
......@@ -28,8 +30,19 @@ class DistanceMap
private:
size_t dim;
std::vector<float> dist;
std::unordered_map<std::string,size_t> index;
struct key_hash
{
size_t operator()(std::tuple<size_t,size_t> const& k) const noexcept
{
size_t i[2];
std::tie(i[0], i[1]) = k;
return boost::hash_range(i, i + 1);
}
};
std::unordered_map<std::tuple<size_t,size_t>,float,key_hash> dist;
};
}
include/cif++/Point.h
......@@ -63,6 +63,15 @@ struct Point
return *this;
}
Point& operator+=(float d)
{
mX += d;
mY += d;
mZ += d;
return *this;
}
Point& operator-=(const Point& rhs)
{
mX -= rhs.mX;
......@@ -72,6 +81,15 @@ struct Point
return *this;
}
Point& operator-=(float d)
{
mX -= d;
mY -= d;
mZ -= d;
return *this;
}
Point& operator*=(float rhs)
{
mX *= rhs;
......
src/BondMap.cpp
......@@ -36,7 +36,7 @@ BondMap::BondMap(const Structure& p)
if (ixb < ixa)
swap(ixa, ixb);
uint32 ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < bond.size());
bond[ix] = true;
......@@ -80,11 +80,11 @@ BondMap::BondMap(const Structure& p)
}
auto c = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == lastSeqID and cif::Key("label_atom_id") == "C");
if (c.size() != 1)
if (c.size() != 1 and VERBOSE > 1)
cerr << "Unexpected number (" << c.size() << ") of atoms with atom ID C in asym_id " << asymId << " with seq id " << lastSeqID << endl;
auto n = db["atom_site"].find(cif::Key("label_asym_id") == asymId and cif::Key("label_seq_id") == seqId and cif::Key("label_atom_id") == "N");
if (n.size() != 1)
if (n.size() != 1 and VERBOSE > 1)
cerr << "Unexpected number (" << n.size() << ") of atoms with atom ID N in asym_id " << asymId << " with seq id " << seqId << endl;
if (not (c.empty() or n.empty()))
......@@ -107,7 +107,7 @@ BondMap::BondMap(const Structure& p)
db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_atom_id") == atomId1) :
db["atom_site"].find(cif::Key("label_asym_id") == asym1 and cif::Key("label_seq_id") == seqId1 and cif::Key("label_atom_id") == atomId1);
if (a.size() != 1)
if (a.size() != 1 and VERBOSE > 1)
cerr << "Unexpected number (" << a.size() << ") of atoms for link with asym_id " << asym1 << " seq_id " << seqId1 << " atom_id " << atomId1 << endl;
auto b =
......@@ -115,7 +115,7 @@ BondMap::BondMap(const Structure& p)
db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_atom_id") == atomId2) :
db["atom_site"].find(cif::Key("label_asym_id") == asym2 and cif::Key("label_seq_id") == seqId2 and cif::Key("label_atom_id") == atomId2);
if (b.size() != 1)
if (b.size() != 1 and VERBOSE > 1)
cerr << "Unexpected number (" << b.size() << ") of atoms for link with asym_id " << asym2 << " seq_id " << seqId2 << " atom_id " << atomId2 << endl;
if (not (a.empty() or b.empty()))
......@@ -186,7 +186,7 @@ BondMap::BondMap(const Structure& p)
if (ixb < ixa)
swap(ixa, ixb);
uint32 ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < bond.size());
bond[ix] = true;
......@@ -201,13 +201,13 @@ BondMap::BondMap(const Structure& p)
bool BondMap::operator()(const Atom& a, const Atom& b) const
{
uint32 ixa = index.at(a.id());
uint32 ixb = index.at(b.id());
size_t ixa = index.at(a.id());
size_t ixb = index.at(b.id());
if (ixb < ixa)
swap(ixa, ixb);
uint32 ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
assert(ix < bond.size());
return bond[ix];
......
src/DistanceMap.cpp
......@@ -3,6 +3,7 @@
#include "cif++/Config.h"
#include <atomic>
#include <mutex>
#include <boost/thread.hpp>
......@@ -44,13 +45,18 @@ vector<clipper::RTop_orth> AlternativeSites(const clipper::Spacegroup& spacegrou
DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegroup, const clipper::Cell& cell)
: dim(0)
{
const float kMaxDistance = 5, kMaxDistanceSQ = kMaxDistance * kMaxDistance;
auto atoms = p.atoms();
dim = atoms.size();
dist = vector<float>(dim * (dim - 1), 0.f);
vector<clipper::Coord_orth> locations(dim);
vector<bool> isWater(dim, false);
// bounding box
Point pMin(numeric_limits<float>::max(), numeric_limits<float>::max(), numeric_limits<float>::max()),
pMax(numeric_limits<float>::min(), numeric_limits<float>::min(), numeric_limits<float>::min());
for (auto& atom: atoms)
{
size_t ix = index.size();
......@@ -58,7 +64,26 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
locations[ix] = atom.location();
isWater[ix] = atom.isWater();
auto p = atom.location();
if (pMin.mX > p.mX)
pMin.mX = p.mX;
if (pMin.mY > p.mY)
pMin.mY = p.mY;
if (pMin.mZ > p.mZ)
pMin.mZ = p.mZ;
if (pMax.mX < p.mX)
pMax.mX = p.mX;
if (pMax.mY < p.mY)
pMax.mY = p.mY;
if (pMax.mZ < p.mZ)
pMax.mZ = p.mZ;
};
pMin -= kMaxDistance; // extend bounding box
pMax += kMaxDistance;
vector<clipper::RTop_orth> rtOrth = AlternativeSites(spacegroup, cell);
......@@ -67,6 +92,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
boost::thread_group t;
size_t N = boost::thread::hardware_concurrency();
atomic<size_t> next(0);
mutex m;
for (size_t i = 0; i < N; ++i)
t.create_thread([&]()
......@@ -77,7 +103,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
if (i >= locations.size())
break;
for (size_t j = i + 1; j < locations.size(); ++j)
{
// if (not (isWater[i] or isWater[j]))
......@@ -90,15 +116,25 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
{
auto p = locations[j].transform(rt);
if (p[0] < pMin.mX or p[1] < pMin.mY or p[2] < pMin.mZ or
p[0] > pMax.mX or p[1] > pMax.mY or p[2] > pMax.mZ)
{
continue;
}
double r2 = (locations[i] - p).lengthsq();
if (minR2 > r2)
minR2 = r2;
}
size_t ix = j + i * dim - i * (i + 1) / 2;
assert(ix < dist.size());
dist[ix] = sqrt(minR2);
if (minR2 < kMaxDistanceSQ)
{
float d = sqrt(minR2);
auto k = make_tuple(i, j);
lock_guard<mutex> lock(m);
dist[k] = d;
}
}
progress.consumed(1);
......@@ -113,8 +149,6 @@ DistanceMap::DistanceMap(const vector<Atom>& atoms)
{
dim = atoms.size();
dist = vector<float>(dim * (dim - 1), 0.f);
for (auto& atom: atoms)
{
size_t ix = index.size();
......@@ -125,9 +159,7 @@ DistanceMap::DistanceMap(const vector<Atom>& atoms)
{
for (size_t j = i + 1; j < dim; ++j)
{
size_t ix = j + i * dim - i * (i + 1) / 2;
assert(ix < dist.size());
dist[ix] = Distance(atoms[i].location(), atoms[j].location());
dist[make_tuple(i, j)] = Distance(atoms[i].location(), atoms[j].location());
}
}
}
......@@ -157,10 +189,16 @@ float DistanceMap::operator()(const Atom& a, const Atom& b) const
if (ixb < ixa)
swap(ixa, ixb);
size_t ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
tuple<size_t,size_t> k{ ixa, ixb };
auto ii = dist.find(k);
float result = 100;
if (ii != dist.end())
result = ii->second;
assert(ix < dist.size());
return dist[ix];
return result;
}
vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
......@@ -185,15 +223,12 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
if (ixb == ixa)
continue;
size_t ix;
if (ixa < ixb)
ix = ixb + ixa * dim - ixa * (ixa + 1) / 2;
else
ix = ixa + ixb * dim - ixb * (ixb + 1) / 2;
auto ii =
ixa < ixb ?
dist.find(make_tuple(ixa, ixb)) :
dist.find(make_tuple(ixb, ixa));
assert(ix < dist.size());
if (dist[ix] != 0 and dist[ix] <= maxDistance)
if (ii != dist.end() and ii->second <= maxDistance)
result.emplace_back(f, i.first);
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment