Skip to content

L
libcifpp
Commit 909a33c0 by maarten
Options

optimisation process started 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@342 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent f66bd7fd
Showing with 514 additions and 464 deletions
include/cif++/Cif++.h
...@@ -236,6 +236,7 @@ namespace detail ...@@ -236,6 +236,7 @@ namespace detail
struct ItemReference struct ItemReference
{ {
const char* mName; const char* mName;
size_t mColumn;
ItemRow* mRow; ItemRow* mRow;
template<typename T> template<typename T>
...@@ -389,14 +390,16 @@ namespace detail ...@@ -389,14 +390,16 @@ namespace detail
return mRow[mColumns[ix]]; return mRow[mColumns[ix]];
} }
getRowResult(const Row& r, C... columns) getRowResult(const Row& r, std::array<size_t, N>&& columns)
: mRow(r), mColumns({{columns...}}) {} : mRow(r), mColumns(std::move(columns))
{
}
template<typename... Ts> template<typename... Ts>
operator std::tuple<Ts...>() const; operator std::tuple<Ts...>() const;
const Row& mRow; const Row& mRow;
std::array<const char*, N> mColumns; std::array<size_t, N> mColumns;
}; };
// we want to be able to tie some variables to a RowResult, for this we use tiewraps // we want to be able to tie some variables to a RowResult, for this we use tiewraps
...@@ -550,30 +553,45 @@ class Row ...@@ -550,30 +553,45 @@ class Row
// TODO: implement real const version? // TODO: implement real const version?
const detail::ItemReference operator[](const char* ItemTag) const const detail::ItemReference operator[](size_t column) const
{
return detail::ItemReference{"<anonymous column>", column, mData};
}
const detail::ItemReference operator[](const char* itemTag) const
{ {
return detail::ItemReference{ItemTag, mData}; size_t column = ColumnForItemTag(itemTag);
return detail::ItemReference{itemTag, column, mData};
} }
detail::ItemReference operator[](const char* ItemTag) detail::ItemReference operator[](const char* itemTag)
{ {
return detail::ItemReference{ItemTag, mData}; size_t column = ColumnForItemTag(itemTag);
return detail::ItemReference{itemTag, column, mData};
} }
const detail::ItemReference operator[](const string& ItemTag) const const detail::ItemReference operator[](const string& itemTag) const
{ {
return detail::ItemReference{ItemTag.c_str(), mData}; size_t column = ColumnForItemTag(itemTag.c_str());
return detail::ItemReference{itemTag.c_str(), column, mData};
} }
detail::ItemReference operator[](const string& ItemTag) detail::ItemReference operator[](const string& itemTag)
{ {
return detail::ItemReference{ItemTag.c_str(), mData}; size_t column = ColumnForItemTag(itemTag.c_str());
return detail::ItemReference{itemTag.c_str(), column, mData};
} }
template<typename... C> template<typename... C>
auto get(C... columns) const -> detail::getRowResult<C...> auto get(C... columns) const -> detail::getRowResult<C...>
{ {
return detail::getRowResult<C...>(*this, columns...); std::array<size_t,sizeof...(C)> cix;
auto c = cix.begin();
for (auto col: { columns... })
*c++ = ColumnForItemTag(col);
return detail::getRowResult<C...>(*this, std::move(cix));
} }
bool operator==(const Row& rhs) const bool operator==(const Row& rhs) const
...@@ -591,9 +609,12 @@ class Row ...@@ -591,9 +609,12 @@ class Row
private: private:
void assign(const string& name, const string& value, bool emplacing); void assign(const string& name, const string& value, bool emplacing);
void assign(size_t column, const string& value, bool emplacing);
void assign(const Item& i, bool emplacing); void assign(const Item& i, bool emplacing);
static void swap(const string& name, ItemRow* a, ItemRow* b); static void swap(size_t column, ItemRow* a, ItemRow* b);
size_t ColumnForItemTag(const char* itemTag) const;
ItemRow* mData; ItemRow* mData;
uint32 mLineNr = 0; uint32 mLineNr = 0;
...@@ -1104,6 +1125,8 @@ class Category ...@@ -1104,6 +1125,8 @@ class Category
iset mandatoryFields() const; iset mandatoryFields() const;
iset keyFields() const; iset keyFields() const;
std::set<size_t> keyFieldsByIndex() const;
void drop(const string& field); void drop(const string& field);
void getTagOrder(vector<string>& tags) const; void getTagOrder(vector<string>& tags) const;
......
include/cif++/CifValidator.h
...@@ -130,7 +130,7 @@ class Validator ...@@ -130,7 +130,7 @@ class Validator
void addCategoryValidator(ValidateCategory&& v); void addCategoryValidator(ValidateCategory&& v);
const ValidateCategory* getValidatorForCategory(std::string category) const; const ValidateCategory* getValidatorForCategory(std::string category) const;
void reportError(const std::string& msg); void reportError(const std::string& msg, bool fatal);
std::string dictName() const { return mName; } std::string dictName() const { return mName; }
void dictName(const std::string& name) { mName = name; } void dictName(const std::string& name) { mName = name; }
......
include/cif++/Structure.h
...@@ -183,18 +183,18 @@ typedef std::vector<Atom> AtomView; ...@@ -183,18 +183,18 @@ typedef std::vector<Atom> AtomView;
class Residue class Residue
{ {
public: public:
Residue() = default; Residue(const Residue& rhs) = delete;
Residue(const Residue& rhs) = default; Residue& operator=(const Residue& rhs) = delete;
Residue& operator=(const Residue& rhs) = default;
Residue(const Structure& structure) Residue(Residue&& rhs);
: mStructure(&structure) {} Residue& operator=(Residue&& rhs);
Residue(const Structure& structure, const std::string& compoundID, Residue(const Structure& structure, const std::string& compoundID,
const std::string& asymID, int seqID = 0, const std::string& asymID, int seqID = 0)
const std::string& altID = "")
: mStructure(&structure), mCompoundID(compoundID) : mStructure(&structure), mCompoundID(compoundID)
, mAsymID(asymID), mAltID(altID), mSeqID(seqID) {} , mAsymID(asymID), mSeqID(seqID) {}
virtual ~Residue() {}
const Compound& compound() const; const Compound& compound() const;
const AtomView& atoms() const; const AtomView& atoms() const;
...@@ -204,7 +204,6 @@ class Residue ...@@ -204,7 +204,6 @@ class Residue
const std::string& compoundID() const { return mCompoundID; } const std::string& compoundID() const { return mCompoundID; }
const std::string& asymID() const { return mAsymID; } const std::string& asymID() const { return mAsymID; }
int seqID() const { return mSeqID; } int seqID() const { return mSeqID; }
const std::string& altID() const { return mAltID; }
int authSeqID() const; int authSeqID() const;
std::string authInsCode() const; std::string authInsCode() const;
...@@ -226,10 +225,14 @@ class Residue ...@@ -226,10 +225,14 @@ class Residue
protected: protected:
Residue() {}
friend class Polymer;
const Structure* mStructure = nullptr; const Structure* mStructure = nullptr;
std::string mCompoundID, mAsymID, mAltID; std::string mCompoundID, mAsymID;
int mSeqID = 0; int mSeqID = 0;
mutable AtomView mAtoms; AtomView mAtoms;
}; };
// -------------------------------------------------------------------- // --------------------------------------------------------------------
...@@ -238,13 +241,16 @@ class Residue ...@@ -238,13 +241,16 @@ class Residue
class Monomer : public Residue class Monomer : public Residue
{ {
public: public:
Monomer(); // Monomer();
Monomer(const Monomer& rhs); Monomer(const Monomer& rhs) = delete;
Monomer& operator=(const Monomer& rhs); Monomer& operator=(const Monomer& rhs) = delete;
Monomer(const Polymer& polymer, uint32 index); Monomer(Monomer&& rhs);
Monomer& operator=(Monomer&& rhs);
// Monomer(const Polymer& polymer, uint32 index);
Monomer(const Polymer& polymer, uint32 index, int seqID, Monomer(const Polymer& polymer, uint32 index, int seqID,
const std::string& compoundID, const std::string& altID); const std::string& compoundID);
// Assuming this is really an amino acid... // Assuming this is really an amino acid...
...@@ -276,62 +282,20 @@ class Monomer : public Residue ...@@ -276,62 +282,20 @@ class Monomer : public Residue
// -------------------------------------------------------------------- // --------------------------------------------------------------------
class Polymer class Polymer : public std::vector<Monomer>
{ {
public: public:
Polymer(const Structure& s, const std::string& asymID); // Polymer(const Structure& s, const std::string& asymID);
Polymer(const Structure& s, const std::string& entityID, const std::string& asymID); Polymer(const Structure& s, const std::string& entityID, const std::string& asymID);
Polymer(const Polymer&) = default;
Polymer& operator=(const Polymer&) = default;
struct iterator : public std::iterator<std::random_access_iterator_tag, Monomer>
{
typedef std::iterator<std::bidirectional_iterator_tag, Monomer> base_type;
typedef base_type::reference reference;
typedef base_type::pointer pointer;
iterator(const Polymer& p, uint32 index); Polymer(const Polymer&) = delete;
iterator(iterator&& rhs); Polymer& operator=(const Polymer&) = delete;
iterator(const iterator& rhs); Polymer(Polymer&& rhs);
iterator& operator=(const iterator& rhs); Polymer& operator=(Polymer&& rhs);
// iterator& operator=(iterator&& rhs);
reference operator*() { return mCurrent; } Monomer& getBySeqID(int seqID);
pointer operator->() { return &mCurrent; } const Monomer& getBySeqID(int seqID) const;
iterator& operator++();
iterator operator++(int)
{
iterator result(*this);
operator++();
return result;
}
iterator& operator--();
iterator operator--(int)
{
iterator result(*this);
operator--();
return result;
}
bool operator==(const iterator& rhs) const { return mPolymer == rhs.mPolymer and mIndex == rhs.mIndex; }
bool operator!=(const iterator& rhs) const { return mPolymer != rhs.mPolymer or mIndex != rhs.mIndex; }
private:
const Polymer* mPolymer;
uint32 mIndex;
Monomer mCurrent;
};
iterator begin() const;
iterator end() const;
size_t size() const { return mPolySeq.size(); }
Monomer operator[](size_t index) const;
Monomer getBySeqID(int seqID) const;
Structure* structure() const { return mStructure; } Structure* structure() const { return mStructure; }
...@@ -344,8 +308,6 @@ class Polymer ...@@ -344,8 +308,6 @@ class Polymer
private: private:
friend struct iterator;
Structure* mStructure; Structure* mStructure;
std::string mEntityID; std::string mEntityID;
std::string mAsymID; std::string mAsymID;
...@@ -398,7 +360,8 @@ class Structure ...@@ -398,7 +360,8 @@ class Structure
const AtomView& atoms() const; const AtomView& atoms() const;
AtomView waters() const; AtomView waters() const;
std::vector<Polymer> polymers() const; const std::vector<Polymer>& polymers() const;
std::vector<Residue> nonPolymers() const; std::vector<Residue> nonPolymers() const;
Atom getAtomById(std::string id) const; Atom getAtomById(std::string id) const;
...@@ -444,7 +407,13 @@ class Structure ...@@ -444,7 +407,13 @@ class Structure
cif::Category& category(const char* name) const; cif::Category& category(const char* name) const;
cif::Datablock& datablock() const; cif::Datablock& datablock() const;
struct StructureImpl* mImpl; void insertCompound(const std::string& compoundID, bool isEntity);
File& mFile;
uint32 mModelNr;
AtomView mAtoms;
std::vector<Polymer> mPolymers;
// std::vector<Residue*> mResidues;
}; };
} }
src/Cif++.cpp
...@@ -224,11 +224,9 @@ const char* ItemReference::c_str() const ...@@ -224,11 +224,9 @@ const char* ItemReference::c_str() const
{ {
// assert(mRow->mCategory); // assert(mRow->mCategory);
auto cix = mRow->mCategory->getColumnIndex(mName);
for (auto iv = mRow->mValues; iv != nullptr; iv = iv->mNext) for (auto iv = mRow->mValues; iv != nullptr; iv = iv->mNext)
{ {
if (iv->mColumnIndex == cix) if (iv->mColumnIndex == mColumn)
{ {
if (iv->mText[0] != '.' or iv->mText[1] != 0) if (iv->mText[0] != '.' or iv->mText[1] != 0)
result = iv->mText; result = iv->mText;
...@@ -247,19 +245,9 @@ const char* ItemReference::c_str(const char* defaultValue) const ...@@ -247,19 +245,9 @@ const char* ItemReference::c_str(const char* defaultValue) const
if (mRow != nullptr and mRow->mCategory != nullptr) if (mRow != nullptr and mRow->mCategory != nullptr)
{ {
// assert(mRow->mCategory);
auto cix = mRow->mCategory->getColumnIndex(mName);
if (cix < mRow->mCategory->mColumns.size())
{
auto iv = mRow->mCategory->mColumns[cix].mValidator;
if (iv != nullptr and not iv->mDefault.empty())
result = iv->mDefault.c_str();
for (auto iv = mRow->mValues; iv != nullptr; iv = iv->mNext) for (auto iv = mRow->mValues; iv != nullptr; iv = iv->mNext)
{ {
if (iv->mColumnIndex == cix) if (iv->mColumnIndex == mColumn)
{ {
// only really non-empty values // only really non-empty values
if (iv->mText[0] != 0 and ((iv->mText[0] != '.' and iv->mText[0] != '?') or iv->mText[1] != 0)) if (iv->mText[0] != 0 and ((iv->mText[0] != '.' and iv->mText[0] != '?') or iv->mText[1] != 0))
...@@ -268,6 +256,12 @@ const char* ItemReference::c_str(const char* defaultValue) const ...@@ -268,6 +256,12 @@ const char* ItemReference::c_str(const char* defaultValue) const
break; break;
} }
} }
if (result == defaultValue and mColumn < mRow->mCategory->mColumns.size()) // not found, perhaps the category has a default defined?
{
auto iv = mRow->mCategory->mColumns[mColumn].mValidator;
if (iv != nullptr and not iv->mDefault.empty())
result = iv->mDefault.c_str();
} }
} }
...@@ -281,7 +275,7 @@ bool ItemReference::empty() const ...@@ -281,7 +275,7 @@ bool ItemReference::empty() const
void ItemReference::swap(ItemReference& b) void ItemReference::swap(ItemReference& b)
{ {
Row::swap(mName, mRow, b.mRow); Row::swap(mColumn, mRow, b.mRow);
} }
} }
...@@ -1196,7 +1190,7 @@ size_t Category::addColumn(const string& name) ...@@ -1196,7 +1190,7 @@ size_t Category::addColumn(const string& name)
{ {
itemValidator = mCatValidator->getValidatorForItem(name); itemValidator = mCatValidator->getValidatorForItem(name);
if (itemValidator == nullptr) if (itemValidator == nullptr)
mValidator->reportError("tag " + name + " not allowed in Category " + mName); mValidator->reportError("tag " + name + " not allowed in Category " + mName, false);
} }
mColumns.push_back({name, itemValidator}); mColumns.push_back({name, itemValidator});
...@@ -1489,7 +1483,8 @@ void Category::getTagOrder(vector<string>& tags) const ...@@ -1489,7 +1483,8 @@ void Category::getTagOrder(vector<string>& tags) const
const detail::ItemReference Category::getFirstItem(const char* itemName) const const detail::ItemReference Category::getFirstItem(const char* itemName) const
{ {
return detail::ItemReference{itemName, mHead}; size_t column = getColumnIndex(itemName);
return detail::ItemReference{itemName, column, mHead};
} }
Category::iterator Category::begin() Category::iterator Category::begin()
...@@ -1518,7 +1513,7 @@ bool Category::isValid() ...@@ -1518,7 +1513,7 @@ bool Category::isValid()
if (mCatValidator == nullptr) if (mCatValidator == nullptr)
{ {
mValidator->reportError("undefined Category " + mName); mValidator->reportError("undefined Category " + mName, false);
return false; return false;
} }
...@@ -1529,7 +1524,7 @@ bool Category::isValid() ...@@ -1529,7 +1524,7 @@ bool Category::isValid()
auto iv = mCatValidator->getValidatorForItem(col.mName); auto iv = mCatValidator->getValidatorForItem(col.mName);
if (iv == nullptr) if (iv == nullptr)
{ {
mValidator->reportError("Field " + col.mName + " is not valid in Category " + mName); mValidator->reportError("Field " + col.mName + " is not valid in Category " + mName, false);
result = false; result = false;
} }
...@@ -1540,7 +1535,7 @@ bool Category::isValid() ...@@ -1540,7 +1535,7 @@ bool Category::isValid()
if (not mandatory.empty()) if (not mandatory.empty())
{ {
mValidator->reportError("In Category " + mName + " the following mandatory fields are missing: " + ba::join(mandatory, ", ")); mValidator->reportError("In Category " + mName + " the following mandatory fields are missing: " + ba::join(mandatory, ", "), false);
result = false; result = false;
} }
...@@ -1569,7 +1564,7 @@ bool Category::isValid() ...@@ -1569,7 +1564,7 @@ bool Category::isValid()
if (iv == nullptr) if (iv == nullptr)
{ {
mValidator->reportError("invalid field " + mColumns[cix].mName + " for Category " + mName); mValidator->reportError("invalid field " + mColumns[cix].mName + " for Category " + mName, false);
result = false; result = false;
continue; continue;
} }
...@@ -1588,7 +1583,7 @@ bool Category::isValid() ...@@ -1588,7 +1583,7 @@ bool Category::isValid()
if (iv != nullptr and iv->mMandatory) if (iv != nullptr and iv->mMandatory)
{ {
mValidator->reportError("missing mandatory field " + mColumns[cix].mName + " for Category " + mName); mValidator->reportError("missing mandatory field " + mColumns[cix].mName + " for Category " + mName, false);
result = false; result = false;
} }
} }
...@@ -1610,7 +1605,7 @@ iset Category::fields() const ...@@ -1610,7 +1605,7 @@ iset Category::fields() const
throw runtime_error("No Validator specified"); throw runtime_error("No Validator specified");
if (mCatValidator == nullptr) if (mCatValidator == nullptr)
mValidator->reportError("undefined Category"); mValidator->reportError("undefined Category", true);
iset result; iset result;
for (auto& iv: mCatValidator->mItemValidators) for (auto& iv: mCatValidator->mItemValidators)
...@@ -1624,7 +1619,7 @@ iset Category::mandatoryFields() const ...@@ -1624,7 +1619,7 @@ iset Category::mandatoryFields() const
throw runtime_error("No Validator specified"); throw runtime_error("No Validator specified");
if (mCatValidator == nullptr) if (mCatValidator == nullptr)
mValidator->reportError("undefined Category"); mValidator->reportError("undefined Category", true);
return mCatValidator->mMandatoryFields; return mCatValidator->mMandatoryFields;
} }
...@@ -1635,11 +1630,26 @@ iset Category::keyFields() const ...@@ -1635,11 +1630,26 @@ iset Category::keyFields() const
throw runtime_error("No Validator specified"); throw runtime_error("No Validator specified");
if (mCatValidator == nullptr) if (mCatValidator == nullptr)
mValidator->reportError("undefined Category"); mValidator->reportError("undefined Category", true);
return iset{ mCatValidator->mKeys.begin(), mCatValidator->mKeys.end() }; return iset{ mCatValidator->mKeys.begin(), mCatValidator->mKeys.end() };
} }
set<size_t> Category::keyFieldsByIndex() const
{
if (mValidator == nullptr)
throw runtime_error("No Validator specified");
if (mCatValidator == nullptr)
mValidator->reportError("undefined Category", true);
set<size_t> result;
for (auto& k: mCatValidator->mKeys)
result.insert(getColumnIndex(k));
return result;
}
auto Category::iterator::operator++() -> iterator& auto Category::iterator::operator++() -> iterator&
{ {
mCurrent = Row(mCurrent.data()->mNext); mCurrent = Row(mCurrent.data()->mNext);
...@@ -1901,12 +1911,24 @@ Row& Row::operator=(const Row& rhs) ...@@ -1901,12 +1911,24 @@ Row& Row::operator=(const Row& rhs)
void Row::assign(const string& name, const string& value, bool emplacing) void Row::assign(const string& name, const string& value, bool emplacing)
{ {
try
{
auto cat = mData->mCategory;
assign(cat->addColumn(name), value, emplacing);
}
catch (const exception& ex)
{
throw_with_nested(logic_error("Could not assigning value '" + value + "' to column " + name));
}
}
void Row::assign(size_t column, const string& value, bool emplacing)
{
if (mData == nullptr) if (mData == nullptr)
throw logic_error("invalid Row, no data assigning value '" + value + "' to " + name); throw logic_error("invalid Row, no data assigning value '" + value + "' to column with index " + to_string(column));
auto cat = mData->mCategory; auto cat = mData->mCategory;
auto cix = cat->addColumn(name); auto& col = cat->mColumns[column];
auto& col = cat->mColumns[cix];
auto& db = cat->mDb; auto& db = cat->mDb;
const char* oldValue = nullptr; const char* oldValue = nullptr;
...@@ -1914,7 +1936,7 @@ void Row::assign(const string& name, const string& value, bool emplacing) ...@@ -1914,7 +1936,7 @@ void Row::assign(const string& name, const string& value, bool emplacing)
{ {
assert(iv != iv->mNext and (iv->mNext == nullptr or iv != iv->mNext->mNext)); assert(iv != iv->mNext and (iv->mNext == nullptr or iv != iv->mNext->mNext));
if (iv->mColumnIndex == cix) if (iv->mColumnIndex == column)
{ {
oldValue = iv->mText; oldValue = iv->mText;
break; break;
...@@ -1936,7 +1958,7 @@ void Row::assign(const string& name, const string& value, bool emplacing) ...@@ -1936,7 +1958,7 @@ void Row::assign(const string& name, const string& value, bool emplacing)
bool reinsert = false; bool reinsert = false;
if (not emplacing and // an update of an Item's value if (not emplacing and // an update of an Item's value
cat->mIndex != nullptr and cat->keyFields().count(name)) cat->mIndex != nullptr and cat->keyFieldsByIndex().count(column))
{ {
reinsert = cat->mIndex->find(mData); reinsert = cat->mIndex->find(mData);
if (reinsert) if (reinsert)
...@@ -1947,7 +1969,7 @@ void Row::assign(const string& name, const string& value, bool emplacing) ...@@ -1947,7 +1969,7 @@ void Row::assign(const string& name, const string& value, bool emplacing)
if (mData->mValues == nullptr) if (mData->mValues == nullptr)
; // nothing to do ; // nothing to do
else if (mData->mValues->mColumnIndex == cix) else if (mData->mValues->mColumnIndex == column)
{ {
auto iv = mData->mValues; auto iv = mData->mValues;
mData->mValues = iv->mNext; mData->mValues = iv->mNext;
...@@ -1958,7 +1980,7 @@ void Row::assign(const string& name, const string& value, bool emplacing) ...@@ -1958,7 +1980,7 @@ void Row::assign(const string& name, const string& value, bool emplacing)
{ {
for (auto iv = mData->mValues; iv->mNext != nullptr; iv = iv->mNext) for (auto iv = mData->mValues; iv->mNext != nullptr; iv = iv->mNext)
{ {
if (iv->mNext->mColumnIndex == cix) if (iv->mNext->mColumnIndex == column)
{ {
auto nv = iv->mNext; auto nv = iv->mNext;
iv->mNext = nv->mNext; iv->mNext = nv->mNext;
...@@ -1972,7 +1994,7 @@ void Row::assign(const string& name, const string& value, bool emplacing) ...@@ -1972,7 +1994,7 @@ void Row::assign(const string& name, const string& value, bool emplacing)
if (not value.empty()) if (not value.empty())
{ {
auto nv = new(value.length()) ItemValue(value.c_str(), cix); auto nv = new(value.length()) ItemValue(value.c_str(), column);
if (mData->mValues == nullptr) if (mData->mValues == nullptr)
mData->mValues = nv; mData->mValues = nv;
...@@ -2012,7 +2034,7 @@ cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag < ...@@ -2012,7 +2034,7 @@ cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag <
} }
} }
void Row::swap(const string& name, ItemRow* a, ItemRow* b) void Row::swap(size_t cix, ItemRow* a, ItemRow* b)
{ {
if (a == nullptr or b == nullptr) if (a == nullptr or b == nullptr)
throw logic_error("invalid Rows in swap"); throw logic_error("invalid Rows in swap");
...@@ -2022,7 +2044,6 @@ void Row::swap(const string& name, ItemRow* a, ItemRow* b) ...@@ -2022,7 +2044,6 @@ void Row::swap(const string& name, ItemRow* a, ItemRow* b)
throw logic_error("Categories not same in swap"); throw logic_error("Categories not same in swap");
auto cat = a->mCategory; auto cat = a->mCategory;
auto cix = cat->addColumn(name);
auto& col = cat->mColumns[cix]; auto& col = cat->mColumns[cix];
auto& db = cat->mDb; auto& db = cat->mDb;
...@@ -2031,7 +2052,7 @@ void Row::swap(const string& name, ItemRow* a, ItemRow* b) ...@@ -2031,7 +2052,7 @@ void Row::swap(const string& name, ItemRow* a, ItemRow* b)
bool reinsert = false; bool reinsert = false;
if (cat->mIndex != nullptr and cat->keyFields().count(name)) if (cat->mIndex != nullptr and cat->keyFieldsByIndex().count(cix))
{ {
reinsert = true; reinsert = true;
cat->mIndex->erase(a); cat->mIndex->erase(a);
...@@ -2146,12 +2167,17 @@ cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag < ...@@ -2146,12 +2167,17 @@ cerr << "fixing linked item " << child->mCategory->mName << '.' << child->mTag <
} }
} }
void Row::assign(const Item& value, bool emplacing) void Row::assign(const Item& value, bool emplacing)
{ {
assign(value.name(), value.value(), emplacing); assign(value.name(), value.value(), emplacing);
} }
size_t Row::ColumnForItemTag(const char* itemTag) const
{
auto cat = mData->mCategory;
return cat->getColumnIndex(itemTag);
}
bool Row::empty() const bool Row::empty() const
{ {
return mData == nullptr or mData->mValues == nullptr; return mData == nullptr or mData->mValues == nullptr;
......
src/CifValidator.cpp
...@@ -252,9 +252,9 @@ ValidateItem* Validator::getValidatorForItem(string tag) const ...@@ -252,9 +252,9 @@ ValidateItem* Validator::getValidatorForItem(string tag) const
return result; return result;
} }
void Validator::reportError(const string& msg) void Validator::reportError(const string& msg, bool fatal)
{ {
if (mStrict) if (mStrict or fatal)
throw ValidationError(msg); throw ValidationError(msg);
else if (VERBOSE) else if (VERBOSE)
cerr << msg << endl; cerr << msg << endl;
......
src/Secondary.cpp
...@@ -153,8 +153,8 @@ const uint32 kHistogramSize = 30; ...@@ -153,8 +153,8 @@ const uint32 kHistogramSize = 30;
struct Res struct Res
{ {
Res(Monomer&& m) Res(const Monomer& m)
: mM(move(m)) : mM(m)
, mType(MapResidue(m.compoundID())) , mType(MapResidue(m.compoundID()))
{ {
for (auto& a: mM.atoms()) for (auto& a: mM.atoms())
...@@ -243,7 +243,7 @@ struct Res ...@@ -243,7 +243,7 @@ struct Res
Res* mNext = nullptr; Res* mNext = nullptr;
Res* mPrev = nullptr; Res* mPrev = nullptr;
Monomer mM; const Monomer& mM;
Point mCAlpha, mC, mN, mO, mH; Point mCAlpha, mC, mN, mO, mH;
...@@ -727,7 +727,7 @@ void CalculateAlphaHelices(vector<Res>& inResidues, bool inPreferPiHelices = tru ...@@ -727,7 +727,7 @@ void CalculateAlphaHelices(vector<Res>& inResidues, bool inPreferPiHelices = tru
void CalculateSecondaryStructure(Structure& s) void CalculateSecondaryStructure(Structure& s)
{ {
auto polymers = s.polymers(); auto& polymers = s.polymers();
size_t nRes = accumulate(polymers.begin(), polymers.end(), size_t nRes = accumulate(polymers.begin(), polymers.end(),
0.0, [](double s, auto& p) { return s + p.size(); }); 0.0, [](double s, auto& p) { return s + p.size(); });
...@@ -737,8 +737,8 @@ void CalculateSecondaryStructure(Structure& s) ...@@ -737,8 +737,8 @@ void CalculateSecondaryStructure(Structure& s)
for (auto& p: polymers) for (auto& p: polymers)
{ {
for (auto m: p) for (auto& m: p)
residues.emplace_back(move(m)); residues.emplace_back(m);
} }
for (size_t i = 0; i + 1 < residues.size(); ++i) for (size_t i = 0; i + 1 < residues.size(); ++i)
...@@ -773,7 +773,7 @@ struct DSSPImpl ...@@ -773,7 +773,7 @@ struct DSSPImpl
DSSPImpl(const Structure& s); DSSPImpl(const Structure& s);
const Structure& mStructure; const Structure& mStructure;
vector<Polymer> mPolymers; const vector<Polymer>& mPolymers;
vector<Res> mResidues; vector<Res> mResidues;
}; };
...@@ -826,8 +826,8 @@ auto start = std::chrono::system_clock::now(); ...@@ -826,8 +826,8 @@ auto start = std::chrono::system_clock::now();
for (auto& p: mPolymers) for (auto& p: mPolymers)
{ {
for (auto m: p) for (auto& m: p)
mResidues.emplace_back(move(m)); mResidues.emplace_back(m);
} }
auto end = std::chrono::system_clock::now(); auto end = std::chrono::system_clock::now();
......
src/Structure.cpp
...@@ -587,6 +587,26 @@ float Atom::radius() const ...@@ -587,6 +587,26 @@ float Atom::radius() const
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// residue // residue
Residue::Residue(Residue&& rhs)
: mStructure(rhs.mStructure), mCompoundID(move(rhs.mCompoundID)), mAsymID(move(rhs.mAsymID))
, mSeqID(rhs.mSeqID), mAtoms(move(rhs.mAtoms))
{
cerr << "move constructor residue" << endl;
rhs.mStructure = nullptr;
}
Residue& Residue::operator=(Residue&& rhs)
{
cerr << "move assignment residue" << endl;
mStructure = rhs.mStructure; rhs.mStructure = nullptr;
mCompoundID = move(rhs.mCompoundID);
mAsymID = move(rhs.mAsymID);
mSeqID = rhs.mSeqID;
mAtoms = move(rhs.mAtoms);
return *this;
}
string Residue::authInsCode() const string Residue::authInsCode() const
{ {
assert(mStructure); assert(mStructure);
...@@ -638,7 +658,9 @@ const AtomView& Residue::atoms() const ...@@ -638,7 +658,9 @@ const AtomView& Residue::atoms() const
if (mAtoms.empty()) if (mAtoms.empty())
{ {
#if 0
//#if 0
for (auto& a: mStructure->atoms()) for (auto& a: mStructure->atoms())
{ {
if (mSeqID > 0 and a.labelSeqId() != mSeqID) if (mSeqID > 0 and a.labelSeqId() != mSeqID)
...@@ -648,27 +670,31 @@ const AtomView& Residue::atoms() const ...@@ -648,27 +670,31 @@ const AtomView& Residue::atoms() const
a.labelCompId() != mCompoundID) a.labelCompId() != mCompoundID)
continue; continue;
if (not mAltID.empty() and a.property<string>("label_alt_id") != mAltID) // if (not mAltID.empty() and a.property<string>("label_alt_id") != mAltID)
continue; // continue;
//
mAtoms.push_back(a); const_cast<AtomView&>(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) if (not mAtoms.empty())
query = move(query) and cif::Key("label_seq_id") == mSeqID; cerr << "loaded atoms for " << labelID() << endl;
if (not mAltID.empty()) //#else
query = move(query) and (cif::Key("label_alt_id") == cif::Empty() or cif::Key("label_alt_id") == mAltID); // auto& atomSites = mStructure->category("atom_site");
//
auto cifAtoms = atomSites.find(move(query)); // auto query = cif::Key("label_asym_id") == mAsymID and cif::Key("label_comp_id") == mCompoundID;
//
for (auto cifAtom: cifAtoms) // if (mSeqID != 0)
mAtoms.push_back(mStructure->getAtomById(cifAtom["id"].as<string>())); // query = move(query) and cif::Key("label_seq_id") == mSeqID;
#endif //
// 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; return mAtoms;
...@@ -688,7 +714,8 @@ const AtomView& Residue::atoms() const ...@@ -688,7 +714,8 @@ const AtomView& Residue::atoms() const
Atom Residue::atomByID(const string& atomID) const Atom Residue::atomByID(const string& atomID) const
{ {
for (auto& a: atoms()) // for (auto& a: atoms())
for (auto& a: mAtoms)
{ {
if (a.labelAtomId() == atomID) if (a.labelAtomId() == atomID)
return a; return a;
...@@ -704,7 +731,8 @@ bool Residue::isEntity() const ...@@ -704,7 +731,8 @@ bool Residue::isEntity() const
auto& db = mStructure->datablock(); auto& db = mStructure->datablock();
auto a1 = db["atom_site"].find(cif::Key("label_asym_id") == mAsymID); auto a1 = db["atom_site"].find(cif::Key("label_asym_id") == mAsymID);
auto a2 = atoms(); // auto a2 = atoms();
auto& a2 = mAtoms;
return a1.size() == a2.size(); return a1.size() == a2.size();
} }
...@@ -740,41 +768,43 @@ string Residue::labelID() const ...@@ -740,41 +768,43 @@ string Residue::labelID() const
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// monomer // monomer
Monomer::Monomer() //Monomer::Monomer(Monomer&& rhs)
: mPolymer(nullptr), mIndex(0) // : Residue(move(rhs)), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex)
{ //{
} //}
Monomer::Monomer(const Monomer& rhs) Monomer::Monomer(const Polymer& polymer, uint32 index, int seqID, const string& compoundID)
: Residue(rhs), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex) : Residue(*polymer.structure(), compoundID, polymer.asymID(), seqID)
, mPolymer(&polymer)
, mIndex(index)
{ {
} }
Monomer& Monomer::operator=(const Monomer& rhs) Monomer::Monomer(Monomer&& rhs)
: Residue(move(rhs)), mPolymer(rhs.mPolymer), mIndex(rhs.mIndex)
{ {
if (this != &rhs) cerr << "move constructor monomer" << endl;
{
Residue::operator=(rhs);
mPolymer = rhs.mPolymer;
mIndex = rhs.mIndex;
}
return *this; // mStructure = rhs.mStructure; rhs.mStructure = nullptr;
// mCompoundID = move(rhs.mCompoundID);
// mAsymID = move(rhs.mAsymID);
// mSeqID = rhs.mSeqID;
// mAtoms = move(rhs.mAtoms);
//
// mPolymer = rhs.mPolymer; rhs.mPolymer = nullptr;
// mIndex = rhs.mIndex;
rhs.mPolymer = nullptr;
} }
Monomer::Monomer(const Polymer& polymer, uint32 index) Monomer& Monomer::operator=(Monomer&& rhs)
: Residue(*polymer.structure())
, mPolymer(&polymer)
, mIndex(index)
{ {
cerr << "move assignment monomer" << endl;
} Residue::operator=(move(rhs));
mPolymer = rhs.mPolymer; rhs.mPolymer = nullptr;
mIndex = rhs.mIndex;
Monomer::Monomer(const Polymer& polymer, uint32 index, int seqID, const string& compoundID, const string& altID) return *this;
: Residue(*polymer.structure(), compoundID, polymer.asymID(), seqID, altID)
, mPolymer(&polymer)
, mIndex(index)
{
} }
float Monomer::phi() const float Monomer::phi() const
...@@ -785,7 +815,7 @@ float Monomer::phi() const ...@@ -785,7 +815,7 @@ float Monomer::phi() const
{ {
if (mIndex > 0) if (mIndex > 0)
{ {
Monomer prev = mPolymer->operator[](mIndex - 1); auto& prev = mPolymer->operator[](mIndex - 1);
if (prev.mSeqID + 1 == mSeqID) if (prev.mSeqID + 1 == mSeqID)
result = DihedralAngle(prev.C().location(), N().location(), CAlpha().location(), C().location()); result = DihedralAngle(prev.C().location(), N().location(), CAlpha().location(), C().location());
} }
...@@ -808,7 +838,7 @@ float Monomer::psi() const ...@@ -808,7 +838,7 @@ float Monomer::psi() const
{ {
if (mIndex + 1 < mPolymer->size()) if (mIndex + 1 < mPolymer->size())
{ {
Monomer next = mPolymer->operator[](mIndex + 1); auto& next = mPolymer->operator[](mIndex + 1);
if (mSeqID + 1 == next.mSeqID) if (mSeqID + 1 == next.mSeqID)
result = DihedralAngle(N().location(), CAlpha().location(), C().location(), next.N().location()); result = DihedralAngle(N().location(), CAlpha().location(), C().location(), next.N().location());
} }
...@@ -830,9 +860,9 @@ float Monomer::alpha() const ...@@ -830,9 +860,9 @@ float Monomer::alpha() const
{ {
if (mIndex > 1 and mIndex + 2 < mPolymer->size()) if (mIndex > 1 and mIndex + 2 < mPolymer->size())
{ {
Monomer prev = mPolymer->operator[](mIndex - 1); auto& prev = mPolymer->operator[](mIndex - 1);
Monomer next = mPolymer->operator[](mIndex + 1); auto& next = mPolymer->operator[](mIndex + 1);
Monomer nextNext = mPolymer->operator[](mIndex + 2); auto& nextNext = mPolymer->operator[](mIndex + 2);
result = DihedralAngle(prev.CAlpha().location(), CAlpha().location(), next.CAlpha().location(), nextNext.CAlpha().location()); result = DihedralAngle(prev.CAlpha().location(), CAlpha().location(), next.CAlpha().location(), nextNext.CAlpha().location());
} }
...@@ -854,8 +884,8 @@ float Monomer::kappa() const ...@@ -854,8 +884,8 @@ float Monomer::kappa() const
{ {
if (mIndex > 2 and mIndex + 2 < mPolymer->size()) if (mIndex > 2 and mIndex + 2 < mPolymer->size())
{ {
Monomer prevPrev = mPolymer->operator[](mIndex - 2); auto& prevPrev = mPolymer->operator[](mIndex - 2);
Monomer nextNext = mPolymer->operator[](mIndex + 2); auto& nextNext = mPolymer->operator[](mIndex + 2);
if (prevPrev.mSeqID + 4 == nextNext.mSeqID) if (prevPrev.mSeqID + 4 == nextNext.mSeqID)
{ {
...@@ -881,7 +911,7 @@ bool Monomer::isCis() const ...@@ -881,7 +911,7 @@ bool Monomer::isCis() const
if (mIndex + 1 < mPolymer->size()) if (mIndex + 1 < mPolymer->size())
{ {
Monomer next = mPolymer->operator[](mIndex + 1); auto& next = mPolymer->operator[](mIndex + 1);
result = Monomer::isCis(*this, next); result = Monomer::isCis(*this, next);
} }
...@@ -935,108 +965,145 @@ bool Monomer::isCis(const mmcif::Monomer& a, const mmcif::Monomer& b) ...@@ -935,108 +965,145 @@ bool Monomer::isCis(const mmcif::Monomer& a, const mmcif::Monomer& b)
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// polymer // polymer
//
//Polymer::iterator::iterator(const Polymer& p, uint32 index)
// : mPolymer(&p), mIndex(index), mCurrent(p, index)
//{
// auto& polySeq = mPolymer->mPolySeq;
//
// if (index < polySeq.size())
// {
// int seqID;
// string asymID, monID;
// cif::tie(asymID, seqID, monID) =
// polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
//
// mCurrent = Monomer(*mPolymer, index, seqID, monID, "");
// }
//}
//
//Monomer Polymer::operator[](size_t index) const
//{
// if (index >= mPolySeq.size())
// throw out_of_range("Invalid index for residue in polymer");
//
// string compoundID;
// int seqID;
//
// auto r = mPolySeq[index];
//
// cif::tie(seqID, compoundID) =
// r.get("seq_id", "mon_id");
//
// return Monomer(const_cast<Polymer&>(*this), index, seqID, compoundID, "");
//}
//
//Polymer::iterator::iterator(const iterator& rhs)
// : mPolymer(rhs.mPolymer), mIndex(rhs.mIndex), mCurrent(rhs.mCurrent)
//{
//}
//
//Polymer::iterator& Polymer::iterator::operator++()
//{
// auto& polySeq = mPolymer->mPolySeq;
//
// if (mIndex < polySeq.size())
// ++mIndex;
//
// if (mIndex < polySeq.size())
// {
// int seqID;
// string asymID, monID;
// cif::tie(asymID, seqID, monID) =
// polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
//
// mCurrent = Monomer(*mPolymer, mIndex, seqID, monID, "");
// }
//
// return *this;
//}
Polymer::iterator::iterator(const Polymer& p, uint32 index) //Polymer::Polymer(const Structure& s, const string& asymID)
: mPolymer(&p), mIndex(index), mCurrent(p, index) // : mStructure(const_cast<Structure*>(&s)), mAsymID(asymID)
{ // , mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID))
auto& polySeq = mPolymer->mPolySeq; //{
// mEntityID = mPolySeq.front()["entity_id"].as<string>();
if (index < polySeq.size()) //
{ //#if DEBUG
int seqID; // for (auto r: mPolySeq)
string asymID, monID; // assert(r["entity_id"] == mEntityID);
cif::tie(asymID, seqID, monID) = //#endif
polySeq[mIndex].get("asym_id", "seq_id", "mon_id"); //
//}
mCurrent = Monomer(*mPolymer, index, seqID, monID, "");
}
}
Monomer Polymer::operator[](size_t index) const Polymer::Polymer(Polymer&& rhs)
: vector<Monomer>(move(rhs))
, mStructure(rhs.mStructure)
, mEntityID(move(rhs.mEntityID)), mAsymID(move(rhs.mAsymID)), mPolySeq(move(rhs.mPolySeq))
{ {
if (index >= mPolySeq.size()) rhs.mStructure = nullptr;
throw out_of_range("Invalid index for residue in polymer");
string compoundID;
int seqID;
auto r = mPolySeq[index];
cif::tie(seqID, compoundID) =
r.get("seq_id", "mon_id");
return Monomer(const_cast<Polymer&>(*this), index, seqID, compoundID, "");
} }
Polymer::iterator::iterator(const iterator& rhs) Polymer& Polymer::operator=(Polymer&& rhs)
: mPolymer(rhs.mPolymer), mIndex(rhs.mIndex), mCurrent(rhs.mCurrent)
{ {
vector<Monomer>::operator=(move(rhs));
mStructure = rhs.mStructure; rhs.mStructure = nullptr;
mEntityID = move(rhs.mEntityID);
mAsymID = move(rhs.mAsymID);
mPolySeq = move(rhs.mPolySeq);
return *this;
} }
Polymer::iterator& Polymer::iterator::operator++() Polymer::Polymer(const Structure& s, const string& entityID, const string& asymID)
: mStructure(const_cast<Structure*>(&s)), mEntityID(entityID), mAsymID(asymID)
, mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID and cif::Key("entity_id") == mEntityID))
{ {
auto& polySeq = mPolymer->mPolySeq; map<uint32,uint32> ix;
if (mIndex < polySeq.size()) for (auto r: mPolySeq)
++mIndex;
if (mIndex < polySeq.size())
{ {
int seqID; int seqID;
string asymID, monID; string compoundID;
cif::tie(asymID, seqID, monID) = cif::tie(seqID, compoundID) = r.get("seq_id", "mon_id");
polySeq[mIndex].get("asym_id", "seq_id", "mon_id");
mCurrent = Monomer(*mPolymer, mIndex, seqID, monID, "");
}
return *this; auto index = size();
}
Polymer::Polymer(const Structure& s, const string& asymID) ix[seqID] = index;
: mStructure(const_cast<Structure*>(&s)), mAsymID(asymID) emplace_back(*this, index, seqID, compoundID);
, mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID)) }
{
mEntityID = mPolySeq.front()["entity_id"].as<string>();
#if DEBUG for (auto atom: s.atoms())
for (auto r: mPolySeq) {
assert(r["entity_id"] == mEntityID); if (atom.labelAsymId() != mAsymID)
#endif continue;
} uint32 index = ix.at(atom.labelSeqId());
Polymer::Polymer(const Structure& s, const string& entityID, const string& asymID) at(index).mAtoms.push_back(atom);
: mStructure(const_cast<Structure*>(&s)), mEntityID(entityID), mAsymID(asymID) }
, mPolySeq(s.category("pdbx_poly_seq_scheme").find(cif::Key("asym_id") == mAsymID and cif::Key("entity_id") == mEntityID))
{
} }
Polymer::iterator Polymer::begin() const string Polymer::chainID() const
{ {
return iterator(*this, 0); return mPolySeq.front()["pdb_strand_id"].as<string>();
} }
Polymer::iterator Polymer::end() const Monomer& Polymer::getBySeqID(int seqID)
{ {
return iterator(*this, mPolySeq.size()); for (auto& m: *this)
} if (m.seqID() == seqID)
return m;
string Polymer::chainID() const throw runtime_error("Monomer with seqID " + to_string(seqID) + " not found in polymer " + mAsymID);
{
return mPolySeq.front()["pdb_strand_id"].as<string>();
} }
Monomer Polymer::getBySeqID(int seqID) const const Monomer& Polymer::getBySeqID(int seqID) const
{ {
for (auto m: *this) for (auto& m: *this)
if (m.seqID() == seqID) if (m.seqID() == seqID)
return m; return m;
if (VERBOSE) throw runtime_error("Monomer with seqID " + to_string(seqID) + " not found in polymer " + mAsymID);
cerr << "Monomer with seqID " << seqID << " not found in polymer " << mAsymID << endl;
return Monomer();
} }
int Polymer::Distance(const Monomer& a, const Monomer& b) const int Polymer::Distance(const Monomer& a, const Monomer& b) const
...@@ -1048,7 +1115,7 @@ int Polymer::Distance(const Monomer& a, const Monomer& b) const ...@@ -1048,7 +1115,7 @@ int Polymer::Distance(const Monomer& a, const Monomer& b) const
int ixa = numeric_limits<int>::max(), ixb = numeric_limits<int>::max(); int ixa = numeric_limits<int>::max(), ixb = numeric_limits<int>::max();
int ix = 0, f = 0; int ix = 0, f = 0;
for (auto m: *this) for (auto& m: *this)
{ {
if (m.seqID() == a.seqID()) if (m.seqID() == a.seqID())
ixa = ix, ++f; ixa = ix, ++f;
...@@ -1169,12 +1236,10 @@ cif::File& File::file() ...@@ -1169,12 +1236,10 @@ cif::File& File::file()
// -------------------------------------------------------------------- // --------------------------------------------------------------------
// Structure // Structure
struct StructureImpl Structure::Structure(File& f, uint32 modelNr)
: mFile(f), mModelNr(modelNr)
{ {
StructureImpl(Structure& s, File& f, uint32 modelNr) auto& db = *mFile.impl().mDb;
: mFile(&f), mModelNr(modelNr)
{
auto& db = *mFile->impl().mDb;
auto& atomCat = db["atom_site"]; auto& atomCat = db["atom_site"];
for (auto& a: atomCat) for (auto& a: atomCat)
...@@ -1185,194 +1250,39 @@ struct StructureImpl ...@@ -1185,194 +1250,39 @@ struct StructureImpl
mAtoms.emplace_back(new AtomImpl(f, a["id"].as<string>(), a)); 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(); }); // sort(mAtoms.begin(), mAtoms.end(), [](auto& a, auto& b) { return a.id() < b.id(); });
}
StructureImpl(const StructureImpl& si)
: mFile(si.mFile), mModelNr(si.mModelNr)
{
mAtoms.reserve(si.mAtoms.size());
for (auto& atom: si.mAtoms)
mAtoms.emplace_back(atom.clone());
}
void removeAtom(Atom& a);
void swapAtoms(Atom& a1, Atom& a2);
void moveAtom(Atom& a, Point p);
void changeResidue(Residue& res, const string& newCompound,
const vector<tuple<string,string>>& remappedAtoms);
void insertCompound(const string& compoundID, bool isEntity); // polymers
File* mFile; auto& polySeqScheme = category("pdbx_poly_seq_scheme");
uint32 mModelNr;
AtomView mAtoms;
};
// --------------------------------------------------------------------
void StructureImpl::insertCompound(const string& compoundID, bool isEntity)
{
auto compound = Compound::create(compoundID);
if (compound == nullptr)
throw runtime_error("Trying to insert unknown compound " + compoundID + " (not found in CCP4 monomers lib)");
cif::Datablock& db = *mFile->impl().mDb;
auto& chemComp = db["chem_comp"];
auto r = chemComp.find(cif::Key("id") == compoundID);
if (r.empty())
{
chemComp.emplace({
{ "id", compoundID },
{ "name", compound->name() },
{ "formula", compound->formula() },
{ "type", compound->type() }
});
}
if (isEntity)
{
auto& pdbxEntityNonpoly = db["pdbx_entity_nonpoly"];
if (pdbxEntityNonpoly.find(cif::Key("comp_id") == compoundID).empty())
{
auto& entity = db["entity"];
string entityID = to_string(entity.size() + 1);
entity.emplace({
{ "id", entityID },
{ "type", "non-polymer" },
{ "pdbx_description", compound->name() },
{ "formula_weight", compound->formulaWeight() }
});
pdbxEntityNonpoly.emplace({
{ "entity_id", entityID },
{ "name", compound->name() },
{ "comp_id", compoundID }
});
}
}
}
void StructureImpl::removeAtom(Atom& a)
{
cif::Datablock& db = *mFile->impl().mDb;
auto& atomSites = db["atom_site"];
for (auto i = atomSites.begin(); i != atomSites.end(); ++i)
{
string id;
cif::tie(id) = i->get("id");
if (id == a.id())
{
atomSites.erase(i);
break;
}
}
mAtoms.erase(remove(mAtoms.begin(), mAtoms.end(), a), mAtoms.end());
}
void StructureImpl::swapAtoms(Atom& a1, Atom& a2)
{
cif::Datablock& db = *mFile->impl().mDb;
auto& atomSites = db["atom_site"];
auto r1 = atomSites.find(cif::Key("id") == a1.id());
auto r2 = atomSites.find(cif::Key("id") == a2.id());
if (r1.size() != 1)
throw runtime_error("Cannot swap atoms since the number of atoms with id " + a1.id() + " is " + to_string(r1.size()));
if (r2.size() != 1)
throw runtime_error("Cannot swap atoms since the number of atoms with id " + a2.id() + " is " + to_string(r2.size()));
auto l1 = r1.front()["label_atom_id"];
auto l2 = r2.front()["label_atom_id"];
std::swap(l1, l2);
auto l3 = r1.front()["auth_atom_id"];
auto l4 = r2.front()["auth_atom_id"];
std::swap(l3, l4);
}
void StructureImpl::moveAtom(Atom& a, Point p)
{
a.location(p);
}
void StructureImpl::changeResidue(Residue& res, const string& newCompound,
const vector<tuple<string,string>>& remappedAtoms)
{
cif::Datablock& db = *mFile->impl().mDb;
string entityID;
// First make sure the compound is already known or insert it.
// And if the residue is an entity, we must make sure it exists
insertCompound(newCompound, res.isEntity());
auto& atomSites = db["atom_site"];
auto atoms = res.atoms();
for (auto& a: remappedAtoms)
{
string a1, a2;
tie(a1, a2) = a;
auto i = find_if(atoms.begin(), atoms.end(), [&](const Atom& a) { return a.labelAtomId() == a1; });
if (i == atoms.end())
{
cerr << "Missing atom for atom ID " << a1 << endl;
continue;
}
auto r = atomSites.find(cif::Key("id") == i->id());
if (r.size() != 1)
continue;
if (a1 != a2)
r.front()["label_atom_id"] = a2;
}
for (auto a: atoms) for (auto& r: polySeqScheme)
{ {
auto r = atomSites.find(cif::Key("id") == a.id()); string asymID, entityID, seqID, monID;
assert(r.size() == 1); cif::tie(asymID, entityID, seqID, monID) =
r.get("asym_id", "entity_id", "seq_id", "mon_id");
if (r.size() != 1)
continue;
r.front()["label_comp_id"] = newCompound; if (mPolymers.empty() or mPolymers.back().asymID() != asymID or mPolymers.back().entityID() != entityID)
if (not entityID.empty()) mPolymers.emplace_back(*this, entityID, asymID);
r.front()["label_entity_id"] = entityID;
} }
} }
// -------------------------------------------------------------------- // Structure(const Structure& s)
// : mFile(s.mFile), mModelNr(s.mModelNr)
Structure::Structure(File& f, uint32 modelNr) // {
: mImpl(new StructureImpl(*this, f, modelNr)) // mAtoms.reserve(si.mAtoms.size());
{ // for (auto& atom: si.mAtoms)
} // mAtoms.emplace_back(atom.clone());
// }
//}
Structure::~Structure() Structure::~Structure()
{ {
delete mImpl;
}
Structure::Structure(const Structure& rhs)
: mImpl(new StructureImpl(*rhs.mImpl))
{
} }
const AtomView& Structure::atoms() const const AtomView& Structure::atoms() const
{ {
return mImpl->mAtoms; return mAtoms;
} }
AtomView Structure::waters() const AtomView Structure::waters() const
...@@ -1395,7 +1305,7 @@ AtomView Structure::waters() const ...@@ -1395,7 +1305,7 @@ AtomView Structure::waters() const
} }
} }
for (auto& a: mImpl->mAtoms) for (auto& a: mAtoms)
{ {
if (a.property<string>("label_entity_id") == waterEntityId) if (a.property<string>("label_entity_id") == waterEntityId)
result.push_back(a); result.push_back(a);
...@@ -1404,23 +1314,24 @@ AtomView Structure::waters() const ...@@ -1404,23 +1314,24 @@ AtomView Structure::waters() const
return result; return result;
} }
vector<Polymer> Structure::polymers() const const vector<Polymer>& Structure::polymers() const
{ {
vector<Polymer> result; // vector<Polymer> result;
//
auto& polySeqScheme = category("pdbx_poly_seq_scheme"); // auto& polySeqScheme = category("pdbx_poly_seq_scheme");
//
for (auto& r: polySeqScheme) // for (auto& r: polySeqScheme)
{ // {
string asymID, entityID, seqID, monID; // string asymID, entityID, seqID, monID;
cif::tie(asymID, entityID, seqID, monID) = // cif::tie(asymID, entityID, seqID, monID) =
r.get("asym_id", "entity_id", "seq_id", "mon_id"); // r.get("asym_id", "entity_id", "seq_id", "mon_id");
//
if (result.empty() or result.back().asymID() != asymID or result.back().entityID() != entityID) // if (result.empty() or result.back().asymID() != asymID or result.back().entityID() != entityID)
result.emplace_back(*this, entityID, asymID); // result.emplace_back(*this, entityID, asymID);
} // }
//
return result; // return result;
return mPolymers;
} }
vector<Residue> Structure::nonPolymers() const vector<Residue> Structure::nonPolymers() const
...@@ -1444,10 +1355,13 @@ vector<Residue> Structure::nonPolymers() const ...@@ -1444,10 +1355,13 @@ vector<Residue> Structure::nonPolymers() const
Atom Structure::getAtomById(string id) const Atom Structure::getAtomById(string id) const
{ {
auto i = lower_bound(mImpl->mAtoms.begin(), mImpl->mAtoms.end(), // auto i = lower_bound(mAtoms.begin(), mAtoms.end(),
id, [](auto& a, auto& b) { return a.id() < b; }); // id, [](auto& a, auto& b) { return a.id() < b; });
auto i = find_if(mAtoms.begin(), mAtoms.end(),
[&id](auto& a) { return a.id() == id; });
if (i == mImpl->mAtoms.end() or i->id() != id) if (i == mAtoms.end() or i->id() != id)
throw out_of_range("Could not find atom with id " + id); throw out_of_range("Could not find atom with id " + id);
return *i; return *i;
...@@ -1455,7 +1369,7 @@ Atom Structure::getAtomById(string id) const ...@@ -1455,7 +1369,7 @@ Atom Structure::getAtomById(string id) const
File& Structure::getFile() const File& Structure::getFile() const
{ {
return *mImpl->mFile; return mFile;
} }
cif::Category& Structure::category(const char* name) const cif::Category& Structure::category(const char* name) const
...@@ -1598,31 +1512,149 @@ tuple<string,int,string> Structure::MapPDBToLabel(const string& asymId, int seqI ...@@ -1598,31 +1512,149 @@ tuple<string,int,string> Structure::MapPDBToLabel(const string& asymId, int seqI
cif::Datablock& Structure::datablock() const cif::Datablock& Structure::datablock() const
{ {
return *mImpl->mFile->impl().mDb; return *mFile.impl().mDb;
} }
// -------------------------------------------------------------------- void Structure::insertCompound(const string& compoundID, bool isEntity)
// actions {
auto compound = Compound::create(compoundID);
if (compound == nullptr)
throw runtime_error("Trying to insert unknown compound " + compoundID + " (not found in CCP4 monomers lib)");
cif::Datablock& db = *mFile.impl().mDb;
auto& chemComp = db["chem_comp"];
auto r = chemComp.find(cif::Key("id") == compoundID);
if (r.empty())
{
chemComp.emplace({
{ "id", compoundID },
{ "name", compound->name() },
{ "formula", compound->formula() },
{ "type", compound->type() }
});
}
if (isEntity)
{
auto& pdbxEntityNonpoly = db["pdbx_entity_nonpoly"];
if (pdbxEntityNonpoly.find(cif::Key("comp_id") == compoundID).empty())
{
auto& entity = db["entity"];
string entityID = to_string(entity.size() + 1);
entity.emplace({
{ "id", entityID },
{ "type", "non-polymer" },
{ "pdbx_description", compound->name() },
{ "formula_weight", compound->formulaWeight() }
});
pdbxEntityNonpoly.emplace({
{ "entity_id", entityID },
{ "name", compound->name() },
{ "comp_id", compoundID }
});
}
}
}
void Structure::removeAtom(Atom& a) void Structure::removeAtom(Atom& a)
{ {
mImpl->removeAtom(a); cif::Datablock& db = *mFile.impl().mDb;
auto& atomSites = db["atom_site"];
for (auto i = atomSites.begin(); i != atomSites.end(); ++i)
{
string id;
cif::tie(id) = i->get("id");
if (id == a.id())
{
atomSites.erase(i);
break;
}
}
mAtoms.erase(remove(mAtoms.begin(), mAtoms.end(), a), mAtoms.end());
} }
void Structure::swapAtoms(Atom& a1, Atom& a2) void Structure::swapAtoms(Atom& a1, Atom& a2)
{ {
mImpl->swapAtoms(a1, a2); cif::Datablock& db = *mFile.impl().mDb;
auto& atomSites = db["atom_site"];
auto r1 = atomSites.find(cif::Key("id") == a1.id());
auto r2 = atomSites.find(cif::Key("id") == a2.id());
if (r1.size() != 1)
throw runtime_error("Cannot swap atoms since the number of atoms with id " + a1.id() + " is " + to_string(r1.size()));
if (r2.size() != 1)
throw runtime_error("Cannot swap atoms since the number of atoms with id " + a2.id() + " is " + to_string(r2.size()));
auto l1 = r1.front()["label_atom_id"];
auto l2 = r2.front()["label_atom_id"];
std::swap(l1, l2);
auto l3 = r1.front()["auth_atom_id"];
auto l4 = r2.front()["auth_atom_id"];
std::swap(l3, l4);
} }
void Structure::moveAtom(Atom& a, Point p) void Structure::moveAtom(Atom& a, Point p)
{ {
mImpl->moveAtom(a, p); a.location(p);
} }
void Structure::changeResidue(Residue& res, const string& newCompound, void Structure::changeResidue(Residue& res, const string& newCompound,
const vector<tuple<string,string>>& remappedAtoms) const vector<tuple<string,string>>& remappedAtoms)
{ {
mImpl->changeResidue(res, newCompound, remappedAtoms); cif::Datablock& db = *mFile.impl().mDb;
string entityID;
// First make sure the compound is already known or insert it.
// And if the residue is an entity, we must make sure it exists
insertCompound(newCompound, res.isEntity());
auto& atomSites = db["atom_site"];
auto atoms = res.atoms();
for (auto& a: remappedAtoms)
{
string a1, a2;
tie(a1, a2) = a;
auto i = find_if(atoms.begin(), atoms.end(), [&](const Atom& a) { return a.labelAtomId() == a1; });
if (i == atoms.end())
{
cerr << "Missing atom for atom ID " << a1 << endl;
continue;
}
auto r = atomSites.find(cif::Key("id") == i->id());
if (r.size() != 1)
continue;
if (a1 != a2)
r.front()["label_atom_id"] = a2;
}
for (auto a: atoms)
{
auto r = atomSites.find(cif::Key("id") == a.id());
assert(r.size() == 1);
if (r.size() != 1)
continue;
r.front()["label_comp_id"] = newCompound;
if (not entityID.empty())
r.front()["label_entity_id"] = entityID;
}
} }
} }
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