Skip to content

L
libcifpp
Commit 51b6c7eb by maarten
Options

meer snelheid 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@345 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent 1cf8e7b7
Showing with 173 additions and 60 deletions
include/cif++/DistanceMap.h
......@@ -46,9 +46,10 @@ class DistanceMap
}
};
std::unordered_map<std::tuple<size_t,size_t>,std::tuple<float,size_t>,key_hash> dist;
Point mD; // needed to move atoms to center
std::vector<clipper::RTop_orth> mRtOrth;
std::vector<std::tuple<float,size_t>> mA;
std::vector<size_t> mIA, mJA;
Point mD; // needed to move atoms to center
std::vector<clipper::RTop_orth> mRtOrth;
};
}
include/cif++/Structure.h
......@@ -389,12 +389,14 @@ class Structure
cif::Datablock& datablock() const;
void insertCompound(const std::string& compoundID, bool isEntity);
void updateAtomIndex();
File& mFile;
uint32 mModelNr;
AtomView mAtoms;
std::vector<size_t> mAtomIndex;
std::list<Polymer> mPolymers;
// std::vector<Residue*> mResidues;
};
}
src/DistanceMap.cpp
......@@ -130,9 +130,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
mD.mY = calculateD(my, cell.b());
mD.mZ = calculateD(mz, cell.c());
if (VERBOSE)
cerr << "Calculated d: " << mD.mX << ", " << mD.mY << " and " << mD.mZ << endl;
if (mD.mX != 0 or mD.mY != 0 or mD.mZ != 0)
{
if (VERBOSE)
......@@ -153,6 +150,8 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
atomic<size_t> next(0);
mutex m;
map<tuple<size_t,size_t>,tuple<float,size_t>> dist;
for (size_t i = 0; i < N; ++i)
t.create_thread([&]()
{
......@@ -185,13 +184,6 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
minR2 = r2;
kb = k;
}
#if defined(DEBUG_VOOR_BART)
if (r2 < 3.5 * 3.5 and not rt.equals(clipper::RTop<>::identity(), 0.1))
cout << "symmetry contact between "
<< atoms[i] << " at " << pi << " and "
<< atoms[j] << " at " << pj << endl;
#endif
}
if (minR2 < kMaxDistanceSQ)
......@@ -209,27 +201,107 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
});
t.join_all();
}
// Store as a sparse CSR compressed matrix
size_t nnz = dist.size();
mA.reserve(nnz);
mIA.reserve(dim + 1);
mJA.reserve(nnz);
size_t lastR = 0;
mIA.push_back(0);
for (auto& di: dist)
{
size_t c, r;
tie(r, c) = di.first;
if (r != lastR) // new row
{
for (size_t ri = lastR; ri < r; ++ri)
mIA.push_back(mA.size());
lastR = r;
}
mA.push_back(di.second);
mJA.push_back(c);
}
for (size_t ri = lastR; ri < dim; ++ri)
mIA.push_back(mA.size());
//DistanceMap::DistanceMap(const Structure& p, const vector<Atom>& atoms)
// : structure(p), dim(0)
//{
// dim = atoms.size();
//// debug code
//
// assert(mIA.size() == dim + 1);
// assert(mA.size() == nnz);
// assert(mJA.size() == nnz);
//
//cerr << "nnz: " << nnz << endl;
//
// for (auto& atom: atoms)
// auto get = [&](size_t i, size_t j)
// {
// size_t ix = index.size();
// index[atom.id()] = ix;
// if (i > j)
// std::swap(i, j);
//
// for (size_t cix = mIA[i]; cix < mIA[i + 1]; ++cix)
// {
// if (mJA[cix] == j)
// return std::get<0>(mA[cix]);
// }
//
// return 100.f;
// };
//
// auto get2 = [&](size_t ixa, size_t ixb)
// {
// if (ixb < ixa)
// std::swap(ixa, ixb);
//
// int32 L = mIA[ixa];
// int32 R = mIA[ixa + 1] - 1;
//
// while (L <= R)
// {
// int32 i = (L + R) / 2;
//
// if (mJA[i] == ixb)
// return std::get<0>(mA[i]);
//
// if (mJA[i] < ixb)
// L = i + 1;
// else
// R = i - 1;
// }
//
// return 100.f;
// };
//
// // test all values
// for (size_t i = 0; i < dim; ++i)
// {
// for (size_t j = i + 1; j < dim; ++j)
// for (size_t j = 0; j < dim; ++j)
// {
// dist[make_tuple(i, j)] = Distance(atoms[i].location(), atoms[j].location());
// float a = get(i, j);
//
// auto ixa = i, ixb = j;
// if (ixb < ixa)
// std::swap(ixa, ixb);
//
// tuple<size_t,size_t> k{ ixa, ixb };
//
// auto ii = dist.find(k);
//
// float b = 100;
//
// if (ii != dist.end())
// b = std::get<0>(ii->second);
//
// assert(a == b);
//
// float c = get2(i, j);
// assert(a == c);
// }
// }
//}
}
float DistanceMap::operator()(const Atom& a, const Atom& b) const
{
......@@ -255,23 +327,28 @@ float DistanceMap::operator()(const Atom& a, const Atom& b) const
if (ixb < ixa)
std::swap(ixa, ixb);
size_t L = mIA[ixa];
size_t R = mIA[ixa + 1] - 1;
tuple<size_t,size_t> k{ ixa, ixb };
while (L <= R)
{
size_t i = (L + R) / 2;
auto ii = dist.find(k);
if (mJA[i] == ixb)
return get<0>(mA[i]);
float result = 100;
if (ii != dist.end())
result = get<0>(ii->second);
if (mJA[i] < ixb)
L = i + 1;
else
R = i - 1;
}
return result;
return 100.f;
}
vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
{
vector<Atom> result;
size_t ixa;
try
{
......@@ -282,27 +359,43 @@ vector<Atom> DistanceMap::near(const Atom& a, float maxDistance) const
throw runtime_error("atom " + a.id() + " not found in distance map");
}
for (auto& di: dist)
set<tuple<size_t,size_t>> bix;
for (size_t bi = mIA[ixa]; bi < mIA[ixa + 1]; ++bi)
{
size_t dixa, dixb;
tie(dixa, dixb) = di.first;
if (dixa != ixa and dixb != ixa)
continue;
float distance;
float d;
size_t rti;
tie(distance, rti) = di.second;
if (distance > maxDistance)
continue;
size_t ixb = dixa;
if (ixb == ixa)
ixb = dixb;
Atom a = structure.getAtomById(rIndex.at(ixb));
tie(d, rti) = mA[bi];
if (d <= maxDistance)
bix.insert(make_tuple(mJA[bi], rti));
}
for (size_t i = 0; i + 1 < dim; ++i)
{
for (size_t j = mIA[i]; j < mIA[i + 1]; ++j)
{
if (mJA[j] != ixa)
continue;
float d;
size_t rti;
tie(d, rti) = mA[j];
if (d > maxDistance)
continue;
bix.insert(make_tuple(i, rti));
}
}
vector<Atom> result;
result.reserve(bix.size());
for (auto& i: bix)
{
size_t ix, rti;
tie(ix, rti) = i;
Atom a = structure.getAtomById(rIndex.at(ix));
result.push_back(a.symmetryCopy(mD, mRtOrth.at(rti)));
}
......
src/Structure.cpp
......@@ -2,6 +2,8 @@
#include "cif++/Structure.h"
#include <numeric>
#include <boost/algorithm/string.hpp>
#include <boost/filesystem/fstream.hpp>
#include <boost/iostreams/filtering_stream.hpp>
......@@ -1259,6 +1261,8 @@ Structure::Structure(File& f, uint32 modelNr)
// sort(mAtoms.begin(), mAtoms.end(), [](auto& a, auto& b) { return a.id() < b.id(); });
updateAtomIndex();
// polymers
auto& polySeqScheme = category("pdbx_poly_seq_scheme");
......@@ -1280,6 +1284,8 @@ Structure::Structure(const Structure& s)
for (auto& atom: s.mAtoms)
mAtoms.emplace_back(atom.clone());
updateAtomIndex();
auto& polySeqScheme = category("pdbx_poly_seq_scheme");
for (auto& r: polySeqScheme)
......@@ -1297,6 +1303,15 @@ Structure::~Structure()
{
}
void Structure::updateAtomIndex()
{
mAtomIndex = vector<size_t>(mAtoms.size());
iota(mAtomIndex.begin(), mAtomIndex.end(), 0);
sort(mAtomIndex.begin(), mAtomIndex.end(), [this](size_t a, size_t b) { return mAtoms[a].id() < mAtoms[b].id(); });
}
const AtomView& Structure::atoms() const
{
return mAtoms;
......@@ -1357,16 +1372,16 @@ vector<Residue> Structure::nonPolymers() const
Atom Structure::getAtomById(string id) const
{
// auto i = lower_bound(mAtoms.begin(), mAtoms.end(),
// id, [](auto& a, auto& b) { return a.id() < b; });
auto i = lower_bound(mAtomIndex.begin(), mAtomIndex.end(),
id, [this](auto& a, auto& b) { return mAtoms[a].id() < b; });
auto i = find_if(mAtoms.begin(), mAtoms.end(),
[&id](auto& a) { return a.id() == id; });
// auto i = find_if(mAtoms.begin(), mAtoms.end(),
// [&id](auto& a) { return a.id() == id; });
if (i == mAtoms.end() or i->id() != id)
if (i == mAtomIndex.end() or mAtoms[*i].id() != id)
throw out_of_range("Could not find atom with id " + id);
return *i;
return mAtoms[*i];
}
File& Structure::getFile() const
......@@ -1580,6 +1595,8 @@ void Structure::removeAtom(Atom& a)
}
mAtoms.erase(remove(mAtoms.begin(), mAtoms.end(), a), mAtoms.end());
updateAtomIndex();
}
void Structure::swapAtoms(Atom& a1, Atom& a2)
......
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