Skip to content

L
libcifpp
Commit c9719f87 by Maarten L. Hekkelman
Options

Merge branch 'develop' into trunk

parents 1c9212c7 123d25f8
Showing with 1151 additions and 934 deletions
CMakeLists.txt
......@@ -25,7 +25,7 @@
cmake_minimum_required(VERSION 3.16)
# set the project name
project(cifpp VERSION 2.0.4 LANGUAGES CXX)
project(cifpp VERSION 3.0.1 LANGUAGES CXX)
list(PREPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake")
......@@ -173,10 +173,6 @@ set(CMAKE_THREAD_PREFER_PTHREAD)
set(THREADS_PREFER_PTHREAD_FLAG)
find_package(Threads)
set(Boost_DETAILED_FAILURE_MSG ON)
if(NOT BUILD_SHARED_LIBS)
set(Boost_USE_STATIC_LIBS ON)
endif()
find_package(Boost 1.70.0 REQUIRED COMPONENTS system iostreams regex program_options)
if(NOT MSVC AND Boost_USE_STATIC_LIBS)
......@@ -186,16 +182,14 @@ endif()
# Create a revision file, containing the current git version info
find_package(Git)
if(GIT_FOUND AND EXISTS "${CMAKE_SOURCE_DIR}/.git")
include(GetGitRevisionDescription)
include(GetGitRevisionDescription)
option(GENERATE_CUSTOM_VERSION "Generate a custom version string" OFF)
if(GIT-NOTFOUND OR HEAD-HASH-NOTFOUND OR NOT GENERATE_CUSTOM_VERSION)
get_git_head_revision(REFSPEC COMMITHASH)
# Generate our own version string
git_describe_working_tree(BUILD_VERSION_STRING --match=build --dirty)
else()
message(WARNING "Git not found, cannot set version info")
SET(BUILD_VERSION_STRING ${PROJECT_VERSION})
endif()
......@@ -256,7 +250,6 @@ set(project_headers
${PROJECT_SOURCE_DIR}/include/cif++/CifUtils.hpp
${PROJECT_SOURCE_DIR}/include/cif++/CifValidator.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Compound.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Matrix.hpp
${PROJECT_SOURCE_DIR}/include/cif++/PDB2Cif.hpp
${PROJECT_SOURCE_DIR}/include/cif++/PDB2CifRemark3.hpp
${PROJECT_SOURCE_DIR}/include/cif++/Point.hpp
......
changelog
Version 3.0.0
- Replaced many strings in the API with string_view for
performance reasons.
- Upgraded mmcif::Structure
- various other small fixes
Version 2.0.4
- Reverted a too strict test when reading cif files.
......
include/cif++/AtomType.hpp
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
*
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
*
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
......@@ -29,150 +29,151 @@
#pragma once
#include <cstdint>
#include <string>
#include <stdexcept>
#include <string>
namespace mmcif
{
enum AtomType : uint8_t
{
Nn = 0, // Unknown
H = 1, // Hydro­gen
He = 2, // He­lium
Li = 3, // Lith­ium
Be = 4, // Beryl­lium
B = 5, // Boron
C = 6, // Carbon
N = 7, // Nitro­gen
O = 8, // Oxy­gen
F = 9, // Fluor­ine
Ne = 10, // Neon
Na = 11, // So­dium
Mg = 12, // Magne­sium
Al = 13, // Alumin­ium
Si = 14, // Sili­con
P = 15, // Phos­phorus
S = 16, // Sulfur
Cl = 17, // Chlor­ine
Ar = 18, // Argon
K = 19, // Potas­sium
Ca = 20, // Cal­cium
Sc = 21, // Scan­dium
Ti = 22, // Tita­nium
V = 23, // Vana­dium
Cr = 24, // Chrom­ium
Mn = 25, // Manga­nese
Fe = 26, // Iron
Co = 27, // Cobalt
Ni = 28, // Nickel
Cu = 29, // Copper
Zn = 30, // Zinc
Ga = 31, // Gallium
Ge = 32, // Germa­nium
As = 33, // Arsenic
Se = 34, // Sele­nium
Br = 35, // Bromine
Kr = 36, // Kryp­ton
Rb = 37, // Rubid­ium
Sr = 38, // Stront­ium
Y = 39, // Yttrium
Zr = 40, // Zirco­nium
Nb = 41, // Nio­bium
Mo = 42, // Molyb­denum
Tc = 43, // Tech­netium
Ru = 44, // Ruthe­nium
Rh = 45, // Rho­dium
Pd = 46, // Pallad­ium
Ag = 47, // Silver
Cd = 48, // Cad­mium
In = 49, // Indium
Sn = 50, // Tin
Sb = 51, // Anti­mony
Te = 52, // Tellurium
I = 53, // Iodine
Xe = 54, // Xenon
Cs = 55, // Cae­sium
Ba = 56, // Ba­rium
La = 57, // Lan­thanum
Hf = 72, // Haf­nium
Ta = 73, // Tanta­lum
W = 74, // Tung­sten
Re = 75, // Rhe­nium
Os = 76, // Os­mium
Ir = 77, // Iridium
Pt = 78, // Plat­inum
Au = 79, // Gold
Hg = 80, // Mer­cury
Tl = 81, // Thallium
Pb = 82, // Lead
Bi = 83, // Bis­muth
Po = 84, // Polo­nium
At = 85, // Asta­tine
Rn = 86, // Radon
Fr = 87, // Fran­cium
Ra = 88, // Ra­dium
Ac = 89, // Actin­ium
Rf = 104, // Ruther­fordium
Db = 105, // Dub­nium
Sg = 106, // Sea­borgium
Bh = 107, // Bohr­ium
Hs = 108, // Has­sium
Mt = 109, // Meit­nerium
Ds = 110, // Darm­stadtium
Rg = 111, // Roent­genium
Cn = 112, // Coper­nicium
Nh = 113, // Nihon­ium
Fl = 114, // Flerov­ium
Mc = 115, // Moscov­ium
Lv = 116, // Liver­morium
Ts = 117, // Tenness­ine
Og = 118, // Oga­nesson
Ce = 58, // Cerium
Pr = 59, // Praseo­dymium
Nd = 60, // Neo­dymium
Pm = 61, // Prome­thium
Sm = 62, // Sama­rium
Eu = 63, // Europ­ium
Gd = 64, // Gadolin­ium
Tb = 65, // Ter­bium
Dy = 66, // Dyspro­sium
Ho = 67, // Hol­mium
Er = 68, // Erbium
Tm = 69, // Thulium
Yb = 70, // Ytter­bium
Lu = 71, // Lute­tium
Th = 90, // Thor­ium
Pa = 91, // Protac­tinium
U = 92, // Ura­nium
Np = 93, // Neptu­nium
Pu = 94, // Pluto­nium
Am = 95, // Ameri­cium
Cm = 96, // Curium
Bk = 97, // Berkel­ium
Cf = 98, // Califor­nium
Es = 99, // Einstei­nium
Fm = 100, // Fer­mium
Md = 101, // Mende­levium
No = 102, // Nobel­ium
Lr = 103, // Lawren­cium
D = 129, // Deuterium
Nn = 0, // Unknown
H = 1, // Hydro­gen
He = 2, // He­lium
Li = 3, // Lith­ium
Be = 4, // Beryl­lium
B = 5, // Boron
C = 6, // Carbon
N = 7, // Nitro­gen
O = 8, // Oxy­gen
F = 9, // Fluor­ine
Ne = 10, // Neon
Na = 11, // So­dium
Mg = 12, // Magne­sium
Al = 13, // Alumin­ium
Si = 14, // Sili­con
P = 15, // Phos­phorus
S = 16, // Sulfur
Cl = 17, // Chlor­ine
Ar = 18, // Argon
K = 19, // Potas­sium
Ca = 20, // Cal­cium
Sc = 21, // Scan­dium
Ti = 22, // Tita­nium
V = 23, // Vana­dium
Cr = 24, // Chrom­ium
Mn = 25, // Manga­nese
Fe = 26, // Iron
Co = 27, // Cobalt
Ni = 28, // Nickel
Cu = 29, // Copper
Zn = 30, // Zinc
Ga = 31, // Gallium
Ge = 32, // Germa­nium
As = 33, // Arsenic
Se = 34, // Sele­nium
Br = 35, // Bromine
Kr = 36, // Kryp­ton
Rb = 37, // Rubid­ium
Sr = 38, // Stront­ium
Y = 39, // Yttrium
Zr = 40, // Zirco­nium
Nb = 41, // Nio­bium
Mo = 42, // Molyb­denum
Tc = 43, // Tech­netium
Ru = 44, // Ruthe­nium
Rh = 45, // Rho­dium
Pd = 46, // Pallad­ium
Ag = 47, // Silver
Cd = 48, // Cad­mium
In = 49, // Indium
Sn = 50, // Tin
Sb = 51, // Anti­mony
Te = 52, // Tellurium
I = 53, // Iodine
Xe = 54, // Xenon
Cs = 55, // Cae­sium
Ba = 56, // Ba­rium
La = 57, // Lan­thanum
Hf = 72, // Haf­nium
Ta = 73, // Tanta­lum
W = 74, // Tung­sten
Re = 75, // Rhe­nium
Os = 76, // Os­mium
Ir = 77, // Iridium
Pt = 78, // Plat­inum
Au = 79, // Gold
Hg = 80, // Mer­cury
Tl = 81, // Thallium
Pb = 82, // Lead
Bi = 83, // Bis­muth
Po = 84, // Polo­nium
At = 85, // Asta­tine
Rn = 86, // Radon
Fr = 87, // Fran­cium
Ra = 88, // Ra­dium
Ac = 89, // Actin­ium
Rf = 104, // Ruther­fordium
Db = 105, // Dub­nium
Sg = 106, // Sea­borgium
Bh = 107, // Bohr­ium
Hs = 108, // Has­sium
Mt = 109, // Meit­nerium
Ds = 110, // Darm­stadtium
Rg = 111, // Roent­genium
Cn = 112, // Coper­nicium
Nh = 113, // Nihon­ium
Fl = 114, // Flerov­ium
Mc = 115, // Moscov­ium
Lv = 116, // Liver­morium
Ts = 117, // Tenness­ine
Og = 118, // Oga­nesson
Ce = 58, // Cerium
Pr = 59, // Praseo­dymium
Nd = 60, // Neo­dymium
Pm = 61, // Prome­thium
Sm = 62, // Sama­rium
Eu = 63, // Europ­ium
Gd = 64, // Gadolin­ium
Tb = 65, // Ter­bium
Dy = 66, // Dyspro­sium
Ho = 67, // Hol­mium
Er = 68, // Erbium
Tm = 69, // Thulium
Yb = 70, // Ytter­bium
Lu = 71, // Lute­tium
Th = 90, // Thor­ium
Pa = 91, // Protac­tinium
U = 92, // Ura­nium
Np = 93, // Neptu­nium
Pu = 94, // Pluto­nium
Am = 95, // Ameri­cium
Cm = 96, // Curium
Bk = 97, // Berkel­ium
Cf = 98, // Califor­nium
Es = 99, // Einstei­nium
Fm = 100, // Fer­mium
Md = 101, // Mende­levium
No = 102, // Nobel­ium
Lr = 103, // Lawren­cium
D = 129, // Deuterium
};
// --------------------------------------------------------------------
// AtomTypeInfo
enum RadiusType {
enum RadiusType
{
eRadiusCalculated,
eRadiusEmpirical,
eRadiusCovalentEmpirical,
......@@ -188,12 +189,12 @@ enum RadiusType {
struct AtomTypeInfo
{
AtomType type;
std::string name;
std::string symbol;
float weight;
bool metal;
float radii[eRadiusTypeCount];
AtomType type;
std::string name;
std::string symbol;
float weight;
bool metal;
float radii[eRadiusTypeCount];
};
extern const AtomTypeInfo kKnownAtoms[];
......@@ -205,25 +206,25 @@ class AtomTypeTraits
{
public:
AtomTypeTraits(AtomType a);
AtomTypeTraits(const std::string& symbol);
AtomType type() const { return mInfo->type; }
std::string name() const { return mInfo->name; }
std::string symbol() const { return mInfo->symbol; }
float weight() const { return mInfo->weight; }
bool isMetal() const { return mInfo->metal; }
static bool isElement(const std::string& symbol);
static bool isMetal(const std::string& symbol);
AtomTypeTraits(const std::string &symbol);
AtomType type() const { return mInfo->type; }
std::string name() const { return mInfo->name; }
std::string symbol() const { return mInfo->symbol; }
float weight() const { return mInfo->weight; }
bool isMetal() const { return mInfo->metal; }
static bool isElement(const std::string &symbol);
static bool isMetal(const std::string &symbol);
float radius(RadiusType type = eRadiusSingleBond) const
{
if (type >= eRadiusTypeCount)
throw std::invalid_argument("invalid radius requested");
return mInfo->radii[type] / 100.f;
}
// data type encapsulating the Waasmaier & Kirfel scattering factors
// in a simplified form (only a and b).
// Added the electrion scattering factors as well
......@@ -231,15 +232,18 @@ class AtomTypeTraits
{
double a[6], b[6];
};
// to get the Cval and Siva values, use this constant as charge:
enum { kWKSFVal = -99 };
const SFData& wksf(int charge = 0) const;
const SFData& elsf() const;
enum
{
kWKSFVal = -99
};
const SFData &wksf(int charge = 0) const;
const SFData &elsf() const;
private:
const struct AtomTypeInfo* mInfo;
const struct AtomTypeInfo *mInfo;
};
}
} // namespace mmcif
include/cif++/BondMap.hpp
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
*
* Copyright (c) 2020 NKI/AVL, Netherlands Cancer Institute
*
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
......@@ -26,9 +26,9 @@
#pragma once
#include <unordered_map>
#include <filesystem>
#include <stdexcept>
#include <unordered_map>
#include "cif++/Structure.hpp"
......@@ -38,39 +38,40 @@ namespace mmcif
class BondMapException : public std::runtime_error
{
public:
BondMapException(const std::string& msg)
: runtime_error(msg) {}
BondMapException(const std::string &msg)
: runtime_error(msg)
{
}
};
class BondMap
{
public:
BondMap(const Structure& p);
BondMap(const BondMap&) = delete;
BondMap& operator=(const BondMap&) = delete;
BondMap(const Structure &p);
BondMap(const BondMap &) = delete;
BondMap &operator=(const BondMap &) = delete;
bool operator()(const Atom& a, const Atom& b) const
bool operator()(const Atom &a, const Atom &b) const
{
return isBonded(index.at(a.id()), index.at(b.id()));
}
bool is1_4(const Atom& a, const Atom& b) const
bool is1_4(const Atom &a, const Atom &b) const
{
uint32_t ixa = index.at(a.id());
uint32_t ixb = index.at(b.id());
return bond_1_4.count(key(ixa, ixb));
}
// links coming from the struct_conn records:
std::vector<std::string> linked(const Atom& a) const;
std::vector<std::string> linked(const Atom &a) const;
// This list of atomID's is comming from either CCD or the CCP4 dictionaries loaded
static std::vector<std::string> atomIDsForCompound(const std::string& compoundID);
private:
static std::vector<std::string> atomIDsForCompound(const std::string &compoundID);
private:
bool isBonded(uint32_t ai, uint32_t bi) const
{
return bond.count(key(ai, bi)) != 0;
......@@ -82,20 +83,19 @@ class BondMap
std::swap(a, b);
return static_cast<uint64_t>(a) | (static_cast<uint64_t>(b) << 32);
}
std::tuple<uint32_t,uint32_t> dekey(uint64_t k) const
std::tuple<uint32_t, uint32_t> dekey(uint64_t k) const
{
return std::make_tuple(
static_cast<uint32_t>(k >> 32),
static_cast<uint32_t>(k)
);
static_cast<uint32_t>(k));
}
uint32_t dim;
std::unordered_map<std::string,uint32_t> index;
std::unordered_map<std::string, uint32_t> index;
std::set<uint64_t> bond, bond_1_4;
std::map<std::string,std::set<std::string>> link;
std::map<std::string, std::set<std::string>> link;
};
}
} // namespace mmcif
include/cif++/Cif++.hpp
......@@ -36,6 +36,10 @@
#include <regex>
#include <set>
#include <sstream>
#include <iomanip>
#include <shared_mutex>
#include <boost/format.hpp>
#include "cif++/CifUtils.hpp"
......@@ -141,14 +145,27 @@ class Item
public:
Item() {}
Item(std::string_view name, char value)
: mName(name)
, mValue({ value })
{
}
template<typename T, std::enable_if_t<std::is_floating_point_v<T>, int> = 0>
Item(std::string_view name, const T& value, const char *fmt)
: mName(name)
, mValue((boost::format(fmt) % value).str())
{
}
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
Item(const std::string &name, const T &value)
Item(const std::string_view name, const T &value)
: mName(name)
, mValue(std::to_string(value))
{
}
Item(const std::string &name, const std::string &value)
Item(const std::string_view name, const std::string_view value)
: mName(name)
, mValue(value)
{
......@@ -221,7 +238,7 @@ class Datablock
using iterator = CategoryList::iterator;
using const_iterator = CategoryList::const_iterator;
Datablock(const std::string &name);
Datablock(const std::string_view name);
~Datablock();
Datablock(const Datablock &) = delete;
......@@ -230,33 +247,31 @@ class Datablock
std::string getName() const { return mName; }
void setName(const std::string &n) { mName = n; }
std::string firstItem(const std::string &tag) const;
iterator begin() { return mCategories.begin(); }
iterator end() { return mCategories.end(); }
const_iterator begin() const { return mCategories.begin(); }
const_iterator end() const { return mCategories.end(); }
Category &operator[](const std::string &name);
Category &operator[](std::string_view name);
std::tuple<iterator, bool> emplace(const std::string &name);
std::tuple<iterator, bool> emplace(std::string_view name);
bool isValid();
void validateLinks() const;
void setValidator(Validator *v);
void setValidator(const Validator *v);
// this one only looks up a Category, returns nullptr if it does not exist
const Category *get(const std::string &name) const;
Category *get(const std::string &name);
const Category *get(std::string_view name) const;
Category *get(std::string_view name);
void getTagOrder(std::vector<std::string> &tags) const;
void write(std::ostream &os, const std::vector<std::string> &order);
void write(std::ostream &os);
// convenience function, add a line to the software category
void add_software(const std::string &name, const std::string &classification,
void add_software(const std::string_view name, const std::string &classification,
const std::string &versionNr, const std::string &versionDate);
friend bool operator==(const Datablock &lhs, const Datablock &rhs);
......@@ -264,9 +279,10 @@ class Datablock
friend std::ostream& operator<<(std::ostream &os, const Datablock &data);
private:
std::list<Category> mCategories;
CategoryList mCategories; // LRU
mutable std::shared_mutex mLock;
std::string mName;
Validator *mValidator;
const Validator *mValidator;
Datablock *mNext;
};
......@@ -350,14 +366,14 @@ namespace detail
private:
friend class ::cif::Row;
ItemReference(const char *name, size_t column, Row &row)
ItemReference(std::string_view name, size_t column, Row &row)
: mName(name)
, mColumn(column)
, mRow(row)
{
}
ItemReference(const char *name, size_t column, const Row &row)
ItemReference(std::string_view name, size_t column, const Row &row)
: mName(name)
, mColumn(column)
, mRow(const_cast<Row &>(row))
......@@ -365,7 +381,7 @@ namespace detail
{
}
const char *mName;
std::string_view mName;
size_t mColumn;
Row &mRow;
bool mConst = false;
......@@ -746,16 +762,16 @@ class Row
return detail::ItemReference(itemTag, column, *this);
}
const detail::ItemReference operator[](const std::string &itemTag) const
const detail::ItemReference operator[](std::string_view itemTag) const
{
size_t column = ColumnForItemTag(itemTag.c_str());
return detail::ItemReference(itemTag.c_str(), column, *this);
size_t column = ColumnForItemTag(itemTag);
return detail::ItemReference(itemTag, column, *this);
}
detail::ItemReference operator[](const std::string &itemTag)
detail::ItemReference operator[](std::string_view itemTag)
{
size_t column = ColumnForItemTag(itemTag.c_str());
return detail::ItemReference(itemTag.c_str(), column, *this);
size_t column = ColumnForItemTag(itemTag);
return detail::ItemReference(itemTag, column, *this);
}
template <typename... Ts, size_t N>
......@@ -776,7 +792,7 @@ class Row
}
void assign(const std::vector<Item> &values);
void assign(const std::string &name, const std::string &value, bool updateLinked);
void assign(std::string_view name, const std::string &value, bool updateLinked);
bool operator==(const Row &rhs) const
{
......@@ -803,7 +819,7 @@ class Row
static void swap(size_t column, ItemRow *a, ItemRow *b);
size_t ColumnForItemTag(const char *itemTag) const;
size_t ColumnForItemTag(std::string_view itemTag) const;
mutable ItemRow *mData;
uint32_t mLineNr = 0;
......@@ -1169,11 +1185,11 @@ struct Empty
struct Key
{
Key(const std::string &itemTag)
explicit Key(const std::string &itemTag)
: mItemTag(itemTag)
{
}
Key(const char *itemTag)
explicit Key(const char *itemTag)
: mItemTag(itemTag)
{
}
......@@ -1816,12 +1832,12 @@ class Category
friend class Row;
friend class detail::ItemReference;
Category(Datablock &db, const std::string &name, Validator *Validator);
Category(Datablock &db, const std::string_view name, const Validator *Validator);
Category(const Category &) = delete;
Category &operator=(const Category &) = delete;
~Category();
const std::string name() const { return mName; }
const std::string &name() const { return mName; }
using iterator = iterator_impl<Row>;
using const_iterator = iterator_impl<const Row>;
......@@ -2064,7 +2080,7 @@ class Category
Datablock &db() { return mDb; }
void setValidator(Validator *v);
void setValidator(const Validator *v);
iset fields() const;
iset mandatoryFields() const;
......@@ -2077,8 +2093,8 @@ class Category
void getTagOrder(std::vector<std::string> &tags) const;
// return index for known column, or the next available column index
size_t getColumnIndex(const std::string &name) const;
bool hasColumn(const std::string &name) const;
size_t getColumnIndex(std::string_view name) const;
bool hasColumn(std::string_view name) const;
const std::string &getColumnName(size_t columnIndex) const;
std::vector<std::string> getColumnNames() const;
......@@ -2119,16 +2135,26 @@ class Category
void write(std::ostream &os, const std::vector<std::string> &order);
void write(std::ostream &os, const std::vector<size_t> &order, bool includeEmptyColumns);
size_t addColumn(const std::string &name);
size_t addColumn(std::string_view name);
struct Linked
{
Category *linked;
const ValidateLink *v;
};
void updateLinks();
Datablock &mDb;
std::string mName;
Validator *mValidator;
const Validator *mValidator;
const ValidateCategory *mCatValidator = nullptr;
std::vector<ItemColumn> mColumns;
ItemRow *mHead;
ItemRow *mTail;
class CatIndex *mIndex;
std::vector<Linked> mParentLinks, mChildLinks;
};
// --------------------------------------------------------------------
......@@ -2162,7 +2188,8 @@ class File
void loadDictionary(); // load the default dictionary, that is mmcifDdl in this case
void loadDictionary(const char *dict); // load one of the compiled in dictionaries
void loadDictionary(std::istream &is); // load dictionary from input stream
void setValidator(const Validator *v);
bool isValid();
void validateLinks() const;
......@@ -2183,8 +2210,8 @@ class File
void append(Datablock *e);
Datablock *get(const std::string &name) const;
Datablock &operator[](const std::string &name);
Datablock *get(std::string_view name) const;
Datablock &operator[](std::string_view name);
struct iterator
{
......@@ -2226,10 +2253,8 @@ class File
void getTagOrder(std::vector<std::string> &tags) const;
private:
void setValidator(Validator *v);
Datablock *mHead;
Validator *mValidator;
const Validator *mValidator;
};
// --------------------------------------------------------------------
......
include/cif++/CifParser.hpp
......@@ -28,8 +28,8 @@
#include "cif++/Cif++.hpp"
#include <stack>
#include <map>
#include <stack>
namespace cif
{
......@@ -39,7 +39,7 @@ namespace cif
class CifParserError : public std::runtime_error
{
public:
CifParserError(uint32_t lineNr, const std::string& message);
CifParserError(uint32_t lineNr, const std::string &message);
};
// --------------------------------------------------------------------
......@@ -48,7 +48,8 @@ extern const uint32_t kMaxLineLength;
extern const uint8_t kCharTraitsTable[128];
enum CharTraitsMask: uint8_t {
enum CharTraitsMask : uint8_t
{
kOrdinaryMask = 1 << 0,
kNonBlankMask = 1 << 1,
kTextLeadMask = 1 << 2,
......@@ -75,13 +76,13 @@ inline bool isTextLead(int ch)
return ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kTextLeadMask) != 0;
}
inline bool isAnyPrint(int ch)
inline bool isAnyPrint(int ch)
{
return ch == '\t' or
(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
return ch == '\t' or
(ch >= 0x20 and ch <= 0x7f and (kCharTraitsTable[ch - 0x20] & kAnyPrintMask) != 0);
}
inline bool isUnquotedString(const char* s)
inline bool isUnquotedString(const char *s)
{
bool result = isOrdinary(*s++);
while (result and *s != 0)
......@@ -94,11 +95,7 @@ inline bool isUnquotedString(const char* s)
// --------------------------------------------------------------------
std::tuple<std::string,std::string> splitTagName(const std::string& tag);
// --------------------------------------------------------------------
using DatablockIndex = std::map<std::string,std::size_t>;
using DatablockIndex = std::map<std::string, std::size_t>;
// --------------------------------------------------------------------
// sac Parser, analogous to SAX Parser (simple api for xml)
......@@ -106,15 +103,15 @@ using DatablockIndex = std::map<std::string,std::size_t>;
class SacParser
{
public:
SacParser(std::istream& is, bool init = true);
SacParser(std::istream &is, bool init = true);
virtual ~SacParser() {}
enum CIFToken
{
eCIFTokenUnknown,
eCIFTokenEOF,
eCIFTokenDATA,
eCIFTokenLOOP,
eCIFTokenGLOBAL,
......@@ -124,7 +121,7 @@ class SacParser
eCIFTokenValue,
};
static const char* kTokenName[];
static const char *kTokenName[];
enum CIFValueType
{
......@@ -137,40 +134,39 @@ class SacParser
eCIFValueUnknown
};
static const char* kValueName[];
static const char *kValueName[];
int getNextChar();
void retract();
void restart();
CIFToken getNextToken();
void match(CIFToken token);
bool parseSingleDatablock(const std::string& datablock);
bool parseSingleDatablock(const std::string &datablock);
DatablockIndex indexDatablocks();
bool parseSingleDatablock(const std::string& datablock, const DatablockIndex &index);
bool parseSingleDatablock(const std::string &datablock, const DatablockIndex &index);
void parseFile();
void parseGlobal();
void parseDataBlock();
virtual void parseSaveFrame();
void parseDictionary();
void error(const std::string& msg);
void error(const std::string &msg);
// production methods, these are pure virtual here
virtual void produceDatablock(const std::string& name) = 0;
virtual void produceCategory(const std::string& name) = 0;
virtual void produceDatablock(const std::string &name) = 0;
virtual void produceCategory(const std::string &name) = 0;
virtual void produceRow() = 0;
virtual void produceItem(const std::string& category, const std::string& item, const std::string& value) = 0;
virtual void produceItem(const std::string &category, const std::string &item, const std::string &value) = 0;
protected:
enum State
{
eStateStart,
......@@ -185,21 +181,21 @@ class SacParser
eStateTextField,
eStateFloat = 100,
eStateInt = 110,
// eStateNumericSuffix = 200,
// eStateNumericSuffix = 200,
eStateValue = 300
};
std::istream& mData;
std::istream &mData;
// Parser state
bool mValidate;
uint32_t mLineNr;
bool mBol;
int mState, mStart;
CIFToken mLookahead;
std::string mTokenValue;
CIFValueType mTokenType;
std::stack<int> mBuffer;
bool mValidate;
uint32_t mLineNr;
bool mBol;
int mState, mStart;
CIFToken mLookahead;
std::string mTokenValue;
CIFValueType mTokenType;
std::stack<int> mBuffer;
};
// --------------------------------------------------------------------
......@@ -207,18 +203,18 @@ class SacParser
class Parser : public SacParser
{
public:
Parser(std::istream& is, File& f, bool init = true);
Parser(std::istream &is, File &f, bool init = true);
virtual void produceDatablock(const std::string& name);
virtual void produceCategory(const std::string& name);
virtual void produceDatablock(const std::string &name);
virtual void produceCategory(const std::string &name);
virtual void produceRow();
virtual void produceItem(const std::string& category, const std::string& item, const std::string& value);
virtual void produceItem(const std::string &category, const std::string &item, const std::string &value);
protected:
File& mFile;
Datablock* mDataBlock;
Datablock::iterator mCat;
Row mRow;
File &mFile;
Datablock *mDataBlock;
Datablock::iterator mCat;
Row mRow;
};
// --------------------------------------------------------------------
......@@ -226,23 +222,21 @@ class Parser : public SacParser
class DictParser : public Parser
{
public:
DictParser(Validator& validator, std::istream& is);
DictParser(Validator &validator, std::istream &is);
~DictParser();
void loadDictionary();
private:
private:
virtual void parseSaveFrame();
bool collectItemTypes();
void linkItems();
Validator& mValidator;
File mFile;
struct DictParserDataImpl* mImpl;
bool mCollectedItemTypes = false;
Validator &mValidator;
File mFile;
struct DictParserDataImpl *mImpl;
bool mCollectedItemTypes = false;
};
}
} // namespace cif
include/cif++/CifUtils.hpp
......@@ -67,8 +67,8 @@ std::string get_version_nr();
// some basic utilities: Since we're using ASCII input only, we define for optimisation
// our own case conversion routines.
bool iequals(const std::string &a, const std::string &b);
int icompare(const std::string &a, const std::string &b);
bool iequals(std::string_view a, std::string_view b);
int icompare(std::string_view a, std::string_view b);
bool iequals(const char *a, const char *b);
int icompare(const char *a, const char *b);
......@@ -100,7 +100,7 @@ inline char tolower(int ch)
// --------------------------------------------------------------------
std::tuple<std::string, std::string> splitTagName(const std::string &tag);
std::tuple<std::string, std::string> splitTagName(std::string_view tag);
// --------------------------------------------------------------------
// generate a cif name, mainly used to generate asym_id's
......
include/cif++/CifValidator.hpp
......@@ -36,18 +36,19 @@
namespace cif
{
struct ValidateCategory;
class ValidatorFactory;
// --------------------------------------------------------------------
class ValidationError : public std::exception
{
public:
ValidationError(const std::string& msg);
ValidationError(const std::string& cat, const std::string& item,
const std::string& msg);
const char* what() const noexcept { return mMsg.c_str(); }
ValidationError(const std::string &msg);
ValidationError(const std::string &cat, const std::string &item,
const std::string &msg);
const char *what() const noexcept { return mMsg.c_str(); }
std::string mMsg;
};
......@@ -55,58 +56,60 @@ class ValidationError : public std::exception
enum class DDL_PrimitiveType
{
Char, UChar, Numb
Char,
UChar,
Numb
};
DDL_PrimitiveType mapToPrimitiveType(const std::string& s);
DDL_PrimitiveType mapToPrimitiveType(std::string_view s);
struct ValidateType
{
std::string mName;
DDL_PrimitiveType mPrimitiveType;
std::string mName;
DDL_PrimitiveType mPrimitiveType;
// std::regex mRx;
boost::regex mRx;
boost::regex mRx;
bool operator<(const ValidateType& rhs) const
bool operator<(const ValidateType &rhs) const
{
return icompare(mName, rhs.mName) < 0;
}
// compare values based on type
// int compare(const std::string& a, const std::string& b) const
// {
// return compare(a.c_str(), b.c_str());
// }
int compare(const char* a, const char* b) const;
// compare values based on type
// int compare(const std::string& a, const std::string& b) const
// {
// return compare(a.c_str(), b.c_str());
// }
int compare(const char *a, const char *b) const;
};
struct ValidateItem
{
std::string mTag;
bool mMandatory;
const ValidateType* mType;
cif::iset mEnums;
std::string mDefault;
bool mDefaultIsNull;
ValidateCategory* mCategory = nullptr;
std::string mTag;
bool mMandatory;
const ValidateType *mType;
cif::iset mEnums;
std::string mDefault;
bool mDefaultIsNull;
ValidateCategory *mCategory = nullptr;
// ItemLinked is used for non-key links
struct ItemLinked
{
ValidateItem* mParent;
std::string mParentItem;
std::string mChildItem;
ValidateItem *mParent;
std::string mParentItem;
std::string mChildItem;
};
std::vector<ItemLinked> mLinked;
bool operator<(const ValidateItem& rhs) const
std::vector<ItemLinked> mLinked;
bool operator<(const ValidateItem &rhs) const
{
return icompare(mTag, rhs.mTag) < 0;
}
bool operator==(const ValidateItem& rhs) const
bool operator==(const ValidateItem &rhs) const
{
return iequals(mTag, rhs.mTag);
}
......@@ -116,22 +119,22 @@ struct ValidateItem
struct ValidateCategory
{
std::string mName;
std::vector<std::string> mKeys;
cif::iset mGroups;
cif::iset mMandatoryFields;
std::set<ValidateItem> mItemValidators;
std::string mName;
std::vector<std::string> mKeys;
cif::iset mGroups;
cif::iset mMandatoryFields;
std::set<ValidateItem> mItemValidators;
bool operator<(const ValidateCategory& rhs) const
bool operator<(const ValidateCategory &rhs) const
{
return icompare(mName, rhs.mName) < 0;
}
void addItemValidator(ValidateItem&& v);
const ValidateItem* getValidatorForItem(std::string tag) const;
const std::set<ValidateItem>& itemValidators() const
void addItemValidator(ValidateItem &&v);
const ValidateItem *getValidatorForItem(std::string_view tag) const;
const std::set<ValidateItem> &itemValidators() const
{
return mItemValidators;
}
......@@ -139,12 +142,12 @@ struct ValidateCategory
struct ValidateLink
{
int mLinkGroupID;
std::string mParentCategory;
std::vector<std::string> mParentKeys;
std::string mChildCategory;
std::vector<std::string> mChildKeys;
std::string mLinkGroupLabel;
int mLinkGroupID;
std::string mParentCategory;
std::vector<std::string> mParentKeys;
std::string mChildCategory;
std::vector<std::string> mChildKeys;
std::string mLinkGroupLabel;
};
// --------------------------------------------------------------------
......@@ -152,47 +155,72 @@ struct ValidateLink
class Validator
{
public:
friend class DictParser;
Validator();
Validator(std::string_view name, std::istream &is);
~Validator();
Validator(const Validator& rhs) = delete;
Validator& operator=(const Validator& rhs) = delete;
Validator(Validator&& rhs);
Validator& operator=(Validator&& rhs);
void addTypeValidator(ValidateType&& v);
const ValidateType* getValidatorForType(std::string typeCode) const;
Validator(const Validator &rhs) = delete;
Validator &operator=(const Validator &rhs) = delete;
void addCategoryValidator(ValidateCategory&& v);
const ValidateCategory* getValidatorForCategory(std::string category) const;
Validator(Validator &&rhs);
Validator &operator=(Validator &&rhs);
void addLinkValidator(ValidateLink&& v);
std::vector<const ValidateLink*> getLinksForParent(const std::string& category) const;
std::vector<const ValidateLink*> getLinksForChild(const std::string& category) const;
friend class DictParser;
friend class ValidatorFactory;
void reportError(const std::string& msg, bool fatal);
std::string dictName() const { return mName; }
void dictName(const std::string& name) { mName = name; }
void addTypeValidator(ValidateType &&v);
const ValidateType *getValidatorForType(std::string_view typeCode) const;
std::string dictVersion() const { return mVersion; }
void dictVersion(const std::string& version) { mVersion = version; }
void addCategoryValidator(ValidateCategory &&v);
const ValidateCategory *getValidatorForCategory(std::string_view category) const;
void addLinkValidator(ValidateLink &&v);
std::vector<const ValidateLink *> getLinksForParent(std::string_view category) const;
std::vector<const ValidateLink *> getLinksForChild(std::string_view category) const;
void reportError(const std::string &msg, bool fatal) const;
std::string dictName() const { return mName; }
void dictName(const std::string &name) { mName = name; }
std::string dictVersion() const { return mVersion; }
void dictVersion(const std::string &version) { mVersion = version; }
private:
// name is fully qualified here:
ValidateItem* getValidatorForItem(std::string name) const;
std::string mName;
std::string mVersion;
bool mStrict = false;
// std::set<uint32_t> mSubCategories;
std::set<ValidateType> mTypeValidators;
std::set<ValidateCategory> mCategoryValidators;
std::vector<ValidateLink> mLinkValidators;
ValidateItem *getValidatorForItem(std::string_view name) const;
std::string mName;
std::string mVersion;
bool mStrict = false;
// std::set<uint32_t> mSubCategories;
std::set<ValidateType> mTypeValidators;
std::set<ValidateCategory> mCategoryValidators;
std::vector<ValidateLink> mLinkValidators;
};
// --------------------------------------------------------------------
class ValidatorFactory
{
public:
static ValidatorFactory &instance()
{
return sInstance;
}
const Validator &operator[](std::string_view dictionary);
private:
static ValidatorFactory sInstance;
ValidatorFactory();
std::mutex mMutex;
std::list<Validator> mValidators;
};
}
} // namespace cif
include/cif++/Compound.hpp
......@@ -104,6 +104,7 @@ class Compound
std::string id() const { return mID; }
std::string name() const { return mName; }
std::string type() const { return mType; }
std::string group() const { return mGroup; }
std::string formula() const { return mFormula; }
float formulaWeight() const { return mFormulaWeight; }
int formalCharge() const { return mFormalCharge; }
......@@ -130,11 +131,12 @@ class Compound
friend class CCP4CompoundFactoryImpl;
Compound(cif::Datablock &db);
Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type);
Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group);
std::string mID;
std::string mName;
std::string mType;
std::string mGroup;
std::string mFormula;
float mFormulaWeight = 0;
int mFormalCharge = 0;
......
include/cif++/Matrix.hpp deleted 100644 → 0
/*-
* SPDX-License-Identifier: BSD-2-Clause
*
* Copyright Maarten L. Hekkelman, Radboud University 2008-2011.
* Copyright (c) 2021 NKI/AVL, Netherlands Cancer Institute
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice, this
* list of conditions and the following disclaimer
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
* DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
* ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
* (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
* LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
* ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
* (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
* SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
// --------------------------------------------------------------------
// uBlas compatible matrix types
#pragma once
#include <iostream>
#include <vector>
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
template <typename T>
class MatrixBase
{
public:
using value_type = T;
virtual ~MatrixBase() {}
virtual uint32_t dim_m() const = 0;
virtual uint32_t dim_n() const = 0;
virtual value_type &operator()(uint32_t i, uint32_t j) { throw std::runtime_error("unimplemented method"); }
virtual value_type operator()(uint32_t i, uint32_t j) const = 0;
MatrixBase &operator*=(const value_type &rhs);
MatrixBase &operator-=(const value_type &rhs);
};
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator*=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) *= rhs;
}
}
return *this;
}
template <typename T>
MatrixBase<T> &MatrixBase<T>::operator-=(const T &rhs)
{
for (uint32_t i = 0; i < dim_m(); ++i)
{
for (uint32_t j = 0; j < dim_n(); ++j)
{
operator()(i, j) -= rhs;
}
}
return *this;
}
template <typename T>
std::ostream &operator<<(std::ostream &lhs, const MatrixBase<T> &rhs)
{
lhs << '[' << rhs.dim_m() << ',' << rhs.dim_n() << ']' << '(';
for (uint32_t i = 0; i < rhs.dim_m(); ++i)
{
lhs << '(';
for (uint32_t j = 0; j < rhs.dim_n(); ++j)
{
if (j > 0)
lhs << ',';
lhs << rhs(i, j);
}
lhs << ')';
}
lhs << ')';
return lhs;
}
template <typename T>
class Matrix : public MatrixBase<T>
{
public:
using value_type = T;
template <typename T2>
Matrix(const MatrixBase<T2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
{
m_data = new value_type[m_m * m_n];
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix()
: m_data(nullptr)
, m_m(0)
, m_n(0)
{
}
Matrix(const Matrix &m)
: m_m(m.m_m)
, m_n(m.m_n)
{
m_data = new value_type[m_m * m_n];
std::copy(m.m_data, m.m_data + (m_m * m_n), m_data);
}
Matrix &operator=(const Matrix &m)
{
value_type *t = new value_type[m.m_m * m.m_n];
std::copy(m.m_data, m.m_data + (m.m_m * m.m_n), t);
delete[] m_data;
m_data = t;
m_m = m.m_m;
m_n = m.m_n;
return *this;
}
Matrix(uint32_t m, uint32_t n, T v = T())
: m_m(m)
, m_n(n)
{
m_data = new value_type[m_m * m_n];
std::fill(m_data, m_data + (m_m * m_n), v);
}
virtual ~Matrix()
{
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_m; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
virtual value_type &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
template <typename Func>
void each(Func f)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
f(m_data[i]);
}
template <typename U>
Matrix &operator/=(U v)
{
for (uint32_t i = 0; i < m_m * m_n; ++i)
m_data[i] /= v;
return *this;
}
private:
value_type *m_data;
uint32_t m_m, m_n;
};
// --------------------------------------------------------------------
template <typename T>
class SymmetricMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
SymmetricMatrix(uint32_t n, T v = T())
: m_owner(true)
, m_n(n)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
m_data = new value_type[N];
std::fill(m_data, m_data + N, v);
}
SymmetricMatrix(const T *data, uint32_t n)
: m_owner(false)
, m_data(const_cast<T *>(data))
, m_n(n)
{
}
virtual ~SymmetricMatrix()
{
if (m_owner)
delete[] m_data;
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
T operator()(uint32_t i, uint32_t j) const;
virtual T &operator()(uint32_t i, uint32_t j);
// erase two rows, add one at the end (for neighbour joining)
void erase_2(uint32_t i, uint32_t j);
template <typename Func>
void each(Func f)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
f(m_data[i]);
}
template <typename U>
SymmetricMatrix &operator/=(U v)
{
uint32_t N = (m_n * (m_n + 1)) / 2;
for (uint32_t i = 0; i < N; ++i)
m_data[i] /= v;
return *this;
}
private:
bool m_owner;
value_type *m_data;
uint32_t m_n;
};
template <typename T>
inline T SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
template <typename T>
inline T &SymmetricMatrix<T>::operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
template <typename T>
void SymmetricMatrix<T>::erase_2(uint32_t di, uint32_t dj)
{
uint32_t s = 0, d = 0;
for (uint32_t i = 0; i < m_n; ++i)
{
for (uint32_t j = 0; j < i; ++j)
{
if (i != di and j != dj and i != dj and j != di)
{
if (s != d)
m_data[d] = m_data[s];
++d;
}
++s;
}
}
--m_n;
}
template <typename T>
class IdentityMatrix : public MatrixBase<T>
{
public:
typedef typename MatrixBase<T>::value_type value_type;
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
virtual uint32_t dim_m() const { return m_n; }
virtual uint32_t dim_n() const { return m_n; }
virtual value_type operator()(uint32_t i, uint32_t j) const
{
value_type result = 0;
if (i == j)
result = 1;
return result;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
for (uint32_t li = 0, rj = 0; li < lhs.dim_m() and rj < rhs.dim_n(); ++li, ++rj)
result(i, j) += lhs(li, j) * rhs(i, rj);
}
}
return result;
}
template <typename T>
Matrix<T> operator*(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs);
result *= rhs;
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, const MatrixBase<T> &rhs)
{
Matrix<T> result(std::min(lhs.dim_m(), rhs.dim_m()), std::min(lhs.dim_n(), rhs.dim_n()));
for (uint32_t i = 0; i < result.dim_m(); ++i)
{
for (uint32_t j = 0; j < result.dim_n(); ++j)
{
result(i, j) = lhs(i, j) - rhs(i, j);
}
}
return result;
}
template <typename T>
Matrix<T> operator-(const MatrixBase<T> &lhs, T rhs)
{
Matrix<T> result(lhs.dim_m(), lhs.dim_n());
result -= rhs;
return result;
}
// template <typename T>
// symmetric_matrix<T> hammingDistance(const MatrixBase<T> &lhs, T rhs);
// template <typename T>
// std::vector<T> sum(const MatrixBase<T> &m);
include/cif++/Secondary.hpp
......@@ -137,6 +137,22 @@ class DSSP
public:
friend class iterator;
ResidueInfo()
: mImpl(nullptr)
{
}
ResidueInfo(const ResidueInfo &rhs)
: mImpl(rhs.mImpl)
{
}
ResidueInfo& operator=(const ResidueInfo &rhs)
{
mImpl = rhs.mImpl;
return *this;
}
explicit operator bool() const { return not empty(); }
bool empty() const { return mImpl == nullptr; }
......
include/cif++/Structure.hpp
......@@ -109,12 +109,12 @@ class Atom
float occupancy() const;
template <typename T>
T property(const std::string &name) const;
T property(const std::string_view name) const;
void property(const std::string &name, const std::string &value);
void property(const std::string_view name, const std::string &value);
template <typename T, std::enable_if_t<std::is_arithmetic_v<T>, int> = 0>
void property(const std::string &name, const T &value)
void property(const std::string_view name, const T &value)
{
property(name, std::to_string(value));
}
......@@ -236,6 +236,8 @@ class Residue
Atom atomByID(const std::string &atomID) const;
const std::string &compoundID() const { return mCompoundID; }
void setCompoundID(const std::string &id) { mCompoundID = id; }
const std::string &asymID() const { return mAsymID; }
int seqID() const { return mSeqID; }
std::string entityID() const;
......@@ -404,7 +406,7 @@ class File : public std::enable_shared_from_this<File>
File(const File &) = delete;
File &operator=(const File &) = delete;
cif::Datablock& createDatablock(const std::string &name);
cif::Datablock& createDatablock(const std::string_view name);
void load(const std::filesystem::path &path);
void save(const std::filesystem::path &path);
......@@ -461,9 +463,24 @@ class Structure
Atom getAtomByLabel(const std::string &atomID, const std::string &asymID,
const std::string &compID, int seqID, const std::string &altID = "");
/// \brief Return the atom closest to point \a p
Atom getAtomByPosition(Point p) const;
/// \brief Return the atom closest to point \a p with atom type \a type in a residue of type \a res_type
Atom getAtomByPositionAndType(Point p, std::string_view type, std::string_view res_type) const;
/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
const Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0) const;
/// \brief Get a residue, if \a seqID is zero, the non-polymers are searched
Residue &getResidue(const std::string &asymID, const std::string &compID, int seqID = 0);
/// \brief Get a the single residue for an asym with id \a asymID
const Residue &getResidue(const std::string &asymID) const;
/// \brief Get a the single residue for an asym with id \a asymID
Residue &getResidue(const std::string &asymID);
// map between auth and label locations
std::tuple<std::string, int, std::string> MapAuthToLabel(const std::string &asymID,
......@@ -488,7 +505,7 @@ class Structure
void removeAtom(Atom &a);
void swapAtoms(Atom &a1, Atom &a2); // swap the labels for these atoms
void moveAtom(Atom &a, Point p); // move atom to a new location
void changeResidue(const Residue &res, const std::string &newCompound,
void changeResidue(Residue &res, const std::string &newCompound,
const std::vector<std::tuple<std::string, std::string>> &remappedAtoms);
/// \brief Create a new non-polymer entity, returns new ID
......@@ -519,15 +536,16 @@ class Structure
void cleanupEmptyCategories();
/// \brief Direct access to underlying data
cif::Category &category(std::string_view name) const;
cif::Datablock &datablock() const;
private:
friend Polymer;
friend Residue;
// friend residue_view;
// friend residue_iterator;
cif::Category &category(const char *name) const;
cif::Datablock &datablock() const;
std::string insertCompound(const std::string &compoundID, bool isEntity);
void loadData();
......
include/cif++/Symmetry.hpp
......@@ -37,6 +37,11 @@ namespace mmcif
// --------------------------------------------------------------------
enum class SpacegroupName
{
full, xHM, Hall
};
struct Spacegroup
{
const char* name;
......@@ -133,6 +138,7 @@ CIFPP_EXPORT extern const std::size_t kSymopNrTableSize;
// --------------------------------------------------------------------
int GetSpacegroupNumber(std::string spacegroup); // alternative for clipper's parsing code
int GetSpacegroupNumber(std::string spacegroup); // alternative for clipper's parsing code, using SpacegroupName::full
int GetSpacegroupNumber(std::string spacegroup, SpacegroupName type); // alternative for clipper's parsing code
}
src/Cif2PDB.cpp
......@@ -249,7 +249,7 @@ size_t WriteContinuedLine(std::ostream& pdbFile, std::string header, int& count,
for (auto& line: lines)
{
ba::to_upper(line);
// ba::to_upper(line);
pdbFile << header;
......@@ -811,7 +811,7 @@ typedef RM<3> RM3;
template<int N>
std::ostream& operator<<(std::ostream& os, RM<N>&& rm)
{
os << "REMARK " << std::setw(3) << std::right << N << " " << rm.mDesc << (rm.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(abs(rm.mWidth)) << std::setprecision(rm.mPrecision);
os << "REMARK " << std::setw(3) << std::right << N << " " << rm.mDesc << (rm.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(std::abs(rm.mWidth)) << std::setprecision(rm.mPrecision);
return os;
}
......@@ -824,7 +824,7 @@ struct SEP
std::ostream& operator<<(std::ostream& os, SEP&& sep)
{
os << sep.mText << (sep.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(abs(sep.mWidth)) << std::setprecision(sep.mPrecision);
os << sep.mText << (sep.mWidth > 0 ? std::left : std::right) << std::fixed << std::setw(std::abs(sep.mWidth)) << std::setprecision(sep.mPrecision);
return os;
}
......
src/CifUtils.cpp
......@@ -123,11 +123,12 @@ const uint8_t kCharToLowerMap[256] =
// --------------------------------------------------------------------
bool iequals(const std::string &a, const std::string &b)
bool iequals(std::string_view a, std::string_view b)
{
bool result = a.length() == b.length();
for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end() and bi != b.end(); ++ai, ++bi)
result = tolower(*ai) == tolower(*bi);
for (auto ai = a.begin(), bi = b.begin(); result and ai != a.end(); ++ai, ++bi)
result = kCharToLowerMap[uint8_t(*ai)] == kCharToLowerMap[uint8_t(*bi)];
// result = tolower(*ai) == tolower(*bi);
return result;
}
......@@ -140,7 +141,7 @@ bool iequals(const char *a, const char *b)
return result and *a == *b;
}
int icompare(const std::string &a, const std::string &b)
int icompare(std::string_view a, std::string_view b)
{
int d = 0;
auto ai = a.begin(), bi = b.begin();
......@@ -193,7 +194,7 @@ std::string toLowerCopy(const std::string &s)
// --------------------------------------------------------------------
std::tuple<std::string, std::string> splitTagName(const std::string &tag)
std::tuple<std::string, std::string> splitTagName(std::string_view tag)
{
if (tag.empty())
throw std::runtime_error("empty tag");
......
src/Compound.cpp
......@@ -128,6 +128,8 @@ Compound::Compound(cif::Datablock &db)
// The name should not contain newline characters since that triggers validation errors later on
ba::replace_all(mName, "\n", "");
mGroup = "non-polymer";
auto &chemCompAtom = db["chem_comp_atom"];
for (auto row : chemCompAtom)
{
......@@ -151,10 +153,11 @@ Compound::Compound(cif::Datablock &db)
}
}
Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type)
Compound::Compound(cif::Datablock &db, const std::string &id, const std::string &name, const std::string &type, const std::string &group)
: mID(id)
, mName(name)
, mType(type)
, mGroup(group)
{
auto &chemCompAtom = db["chem_comp_atom"];
for (auto row : chemCompAtom)
......@@ -411,7 +414,7 @@ CompoundFactoryImpl::CompoundFactoryImpl(const std::filesystem::path &file, std:
auto &db = cifFile["comp_" + id];
mCompounds.push_back(new Compound(db, id, name, type));
mCompounds.push_back(new Compound(db, id, name, type, group));
}
}
else
......@@ -617,7 +620,7 @@ Compound *CCP4CompoundFactoryImpl::create(const std::string &id)
else
type = "non-polymer";
mCompounds.push_back(new Compound(db, id, name, type));
mCompounds.push_back(new Compound(db, id, name, type, group));
result = mCompounds.back();
}
}
......
src/Point.cpp
......@@ -28,12 +28,249 @@
#include <valarray>
#include "cif++/Point.hpp"
#include "cif++/Matrix.hpp"
namespace mmcif
{
// --------------------------------------------------------------------
// We're using expression templates here
template <typename M>
class MatrixExpression
{
public:
uint32_t dim_m() const { return static_cast<const M&>(*this).dim_m(); }
uint32_t dim_n() const { return static_cast<const M&>(*this).dim_n(); }
double &operator()(uint32_t i, uint32_t j)
{
return static_cast<M&>(*this).operator()(i, j);
}
double operator()(uint32_t i, uint32_t j) const
{
return static_cast<const M&>(*this).operator()(i, j);
}
};
// --------------------------------------------------------------------
// matrix is m x n, addressing i,j is 0 <= i < m and 0 <= j < n
// element m i,j is mapped to [i * n + j] and thus storage is row major
class Matrix : public MatrixExpression<Matrix>
{
public:
template <typename M2>
Matrix(const MatrixExpression<M2> &m)
: m_m(m.dim_m())
, m_n(m.dim_n())
, m_data(m_m * m_n)
{
for (uint32_t i = 0; i < m_m; ++i)
{
for (uint32_t j = 0; j < m_n; ++j)
operator()(i, j) = m(i, j);
}
}
Matrix(size_t m, size_t n, double v = 0)
: m_m(m)
, m_n(n)
, m_data(m_m * m_n)
{
std::fill(m_data.begin(), m_data.end(), v);
}
Matrix() = default;
Matrix(Matrix &&m) = default;
Matrix(const Matrix &m) = default;
Matrix &operator=(Matrix &&m) = default;
Matrix &operator=(const Matrix &m) = default;
uint32_t dim_m() const { return m_m; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
double &operator()(uint32_t i, uint32_t j)
{
assert(i < m_m);
assert(j < m_n);
return m_data[i * m_n + j];
}
private:
uint32_t m_m = 0, m_n = 0;
std::vector<double> m_data;
};
// --------------------------------------------------------------------
class SymmetricMatrix : public MatrixExpression<SymmetricMatrix>
{
public:
SymmetricMatrix(uint32_t n, double v = 0)
: m_n(n)
, m_data((m_n * (m_n + 1)) / 2)
{
std::fill(m_data.begin(), m_data.end(), v);
}
SymmetricMatrix() = default;
SymmetricMatrix(SymmetricMatrix &&m) = default;
SymmetricMatrix(const SymmetricMatrix &m) = default;
SymmetricMatrix &operator=(SymmetricMatrix &&m) = default;
SymmetricMatrix &operator=(const SymmetricMatrix &m) = default;
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i < j
? m_data[(j * (j + 1)) / 2 + i]
: m_data[(i * (i + 1)) / 2 + j];
}
double &operator()(uint32_t i, uint32_t j)
{
if (i > j)
std::swap(i, j);
assert(j < m_n);
return m_data[(j * (j + 1)) / 2 + i];
}
private:
uint32_t m_n;
std::vector<double> m_data;
};
class IdentityMatrix : public MatrixExpression<IdentityMatrix>
{
public:
IdentityMatrix(uint32_t n)
: m_n(n)
{
}
uint32_t dim_m() const { return m_n; }
uint32_t dim_n() const { return m_n; }
double operator()(uint32_t i, uint32_t j) const
{
return i == j ? 1 : 0;
}
private:
uint32_t m_n;
};
// --------------------------------------------------------------------
// matrix functions, implemented as expression templates
template<typename M1, typename M2>
class MatrixSubtraction : public MatrixExpression<MatrixSubtraction<M1, M2>>
{
public:
MatrixSubtraction(const M1 &m1, const M2 &m2)
: m_m1(m1), m_m2(m2)
{
assert(m_m1.dim_m() == m_m2.dim_m());
assert(m_m1.dim_n() == m_m2.dim_n());
}
uint32_t dim_m() const { return m_m1.dim_m(); }
uint32_t dim_n() const { return m_m1.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m1(i, j) - m_m2(i, j);
}
private:
const M1 &m_m1;
const M2 &m_m2;
};
template<typename M1, typename M2>
MatrixSubtraction<M1, M2> operator-(const MatrixExpression<M1> &m1, const MatrixExpression<M2> &m2)
{
return MatrixSubtraction(*static_cast<const M1*>(&m1), *static_cast<const M2*>(&m2));
}
template<typename M>
class MatrixMultiplication : public MatrixExpression<MatrixMultiplication<M>>
{
public:
MatrixMultiplication(const M &m, double v)
: m_m(m), m_v(v)
{
}
uint32_t dim_m() const { return m_m.dim_m(); }
uint32_t dim_n() const { return m_m.dim_n(); }
double operator()(uint32_t i, uint32_t j) const
{
return m_m(i, j) * m_v;
}
private:
const M &m_m;
double m_v;
};
template<typename M>
MatrixMultiplication<M> operator*(const MatrixExpression<M> &m, double v)
{
return MatrixMultiplication(*static_cast<const M*>(&m), v);
}
// --------------------------------------------------------------------
template<class M1>
Matrix Cofactors(const M1& m)
{
Matrix cf(m.dim_m(), m.dim_m());
const size_t ixs[4][3] =
{
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
for (size_t x = 0; x < 4; ++x)
{
const size_t* ix = ixs[x];
for (size_t y = 0; y < 4; ++y)
{
const size_t* iy = ixs[y];
cf(x, y) =
m(ix[0], iy[0]) * m(ix[1], iy[1]) * m(ix[2], iy[2]) +
m(ix[0], iy[1]) * m(ix[1], iy[2]) * m(ix[2], iy[0]) +
m(ix[0], iy[2]) * m(ix[1], iy[0]) * m(ix[2], iy[1]) -
m(ix[0], iy[2]) * m(ix[1], iy[1]) * m(ix[2], iy[0]) -
m(ix[0], iy[1]) * m(ix[1], iy[0]) * m(ix[2], iy[2]) -
m(ix[0], iy[0]) * m(ix[1], iy[2]) * m(ix[2], iy[1]);
if ((x + y) % 2 == 1)
cf(x, y) *= -1;
}
}
return cf;
}
// --------------------------------------------------------------------
Quaternion Normalize(Quaternion q)
{
......@@ -102,14 +339,14 @@ Point CenterPoints(std::vector<Point>& Points)
return t;
}
Point Centroid(std::vector<Point>& Points)
Point Centroid(const std::vector<Point>& pts)
{
Point result;
for (Point& pt : Points)
for (auto &pt : pts)
result += pt;
result /= static_cast<float>(Points.size());
result /= static_cast<float>(pts.size());
return result;
}
......@@ -175,7 +412,7 @@ double LargestDepressedQuarticSolution(double a, double b, double c)
Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& pb)
{
// First calculate M, a 3x3 Matrix containing the sums of products of the coordinates of A and B
Matrix<double> M(3, 3, 0);
Matrix M(3, 3, 0);
for (uint32_t i = 0; i < pa.size(); ++i)
{
......@@ -188,7 +425,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
}
// Now calculate N, a symmetric 4x4 Matrix
SymmetricMatrix<double> N(4);
SymmetricMatrix N(4);
N(0, 0) = M(0, 0) + M(1, 1) + M(2, 2);
N(0, 1) = M(1, 2) - M(2, 1);
......@@ -225,6 +462,7 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
M(1, 2) * M(2, 0) * M(0, 1) +
M(2, 1) * M(1, 0) * M(0, 2));
// E is the determinant of N:
double E =
(N(0,0) * N(1,1) - N(0,1) * N(0,1)) * (N(2,2) * N(3,3) - N(2,3) * N(2,3)) +
(N(0,1) * N(0,2) - N(0,0) * N(2,1)) * (N(2,1) * N(3,3) - N(2,3) * N(1,3)) +
......@@ -234,47 +472,22 @@ Quaternion AlignPoints(const std::vector<Point>& pa, const std::vector<Point>& p
(N(0,2) * N(1,3) - N(2,1) * N(0,3)) * (N(0,2) * N(1,3) - N(2,1) * N(0,3));
// solve quartic
double lm = LargestDepressedQuarticSolution(C, D, E);
double lambda = LargestDepressedQuarticSolution(C, D, E);
// calculate t = (N - λI)
Matrix<double> li = IdentityMatrix<double>(4) * lm;
Matrix<double> t = N - li;
Matrix t = N - IdentityMatrix(4) * lambda;
// calculate a Matrix of cofactors for t
Matrix<double> cf(4, 4);
Matrix cf = Cofactors(t);
const uint32_t ixs[4][3] =
int maxR = 0;
for (int r = 1; r < 4; ++r)
{
{ 1, 2, 3 },
{ 0, 2, 3 },
{ 0, 1, 3 },
{ 0, 1, 2 }
};
uint32_t maxR = 0;
for (uint32_t r = 0; r < 4; ++r)
{
const uint32_t* ir = ixs[r];
for (uint32_t c = 0; c < 4; ++c)
{
const uint32_t* ic = ixs[c];
cf(r, c) =
t(ir[0], ic[0]) * t(ir[1], ic[1]) * t(ir[2], ic[2]) +
t(ir[0], ic[1]) * t(ir[1], ic[2]) * t(ir[2], ic[0]) +
t(ir[0], ic[2]) * t(ir[1], ic[0]) * t(ir[2], ic[1]) -
t(ir[0], ic[2]) * t(ir[1], ic[1]) * t(ir[2], ic[0]) -
t(ir[0], ic[1]) * t(ir[1], ic[0]) * t(ir[2], ic[2]) -
t(ir[0], ic[0]) * t(ir[1], ic[2]) * t(ir[2], ic[1]);
}
if (r > maxR and cf(r, 0) > cf(maxR, 0))
if (std::abs(cf(r, 0)) > std::abs(cf(maxR, 0)))
maxR = r;
}
// NOTE the negation of the y here, why? Maybe I swapped r/c above?
Quaternion q(cf(maxR, 0), cf(maxR, 1), -cf(maxR, 2), cf(maxR, 3));
Quaternion q(cf(maxR, 0), cf(maxR, 1), cf(maxR, 2), cf(maxR, 3));
q = Normalize(q);
return q;
......
src/Secondary.cpp
......@@ -897,7 +897,6 @@ void CalculateBetaSheets(std::vector<Res>& inResidues, DSSP_Statistics& stats)
}
// --------------------------------------------------------------------
// TODO: improve alpha helix calculation by better recognizing pi-helices
void CalculateAlphaHelices(std::vector<Res>& inResidues, DSSP_Statistics& stats, bool inPreferPiHelices = true)
{
......
src/Structure.cpp
......@@ -216,9 +216,9 @@ struct AtomImpl
, mLocation(i.mLocation)
, mRefcount(1)
, mRow(i.mRow)
, mCachedRefs(i.mCachedRefs)
, mCompound(i.mCompound)
, mRadius(i.mRadius)
, mCachedProperties(i.mCachedProperties)
, mSymmetryCopy(i.mSymmetryCopy)
, mClone(true)
// , mRTop(i.mRTop), mD(i.mD)
......@@ -270,9 +270,9 @@ struct AtomImpl
, mLocation(loc)
, mRefcount(1)
, mRow(impl.mRow)
, mCachedRefs(impl.mCachedRefs)
, mCompound(impl.mCompound)
, mRadius(impl.mRadius)
, mCachedProperties(impl.mCachedProperties)
, mSymmetryCopy(true)
, mSymmetryOperator(sym_op)
{
......@@ -317,13 +317,15 @@ struct AtomImpl
auto cat = mDb.get("atom_site_anisotrop");
if (cat)
{
auto r = cat->find1(cif::Key("id") == mID);
if (not r.empty())
try
{
result = true;
auto r = cat->find1(cif::Key("id") == mID);
cif::tie(anisou[0], anisou[1], anisou[2], anisou[3], anisou[4], anisou[5]) =
r.get("U[1][1]", "U[1][2]", "U[1][3]", "U[2][2]", "U[2][3]", "U[3][3]");
result = true;
}
catch(const std::exception& e)
{
}
}
......@@ -338,9 +340,9 @@ struct AtomImpl
if (not mClone)
{
mRow["Cartn_x"] = p.getX();
mRow["Cartn_y"] = p.getY();
mRow["Cartn_z"] = p.getZ();
property("Cartn_x", std::to_string(p.getX()));
property("Cartn_y", std::to_string(p.getY()));
property("Cartn_z", std::to_string(p.getZ()));
}
// boost::format kPosFmt("%.3f");
......@@ -382,26 +384,31 @@ struct AtomImpl
return mRadius;
}
const std::string &property(const std::string &name) const
const std::string property(const std::string_view name) const
{
static std::string kEmptyString;
auto i = mCachedProperties.find(name);
if (i == mCachedProperties.end())
for (auto &&[tag, ref] : mCachedRefs)
{
auto v = mRow[name];
if (v.empty())
return kEmptyString;
return mCachedProperties[name] = v.as<std::string>();
if (tag == name)
return ref.as<std::string>();
}
else
return i->second;
mCachedRefs.emplace_back(name, mRow[name]);
return std::get<1>(mCachedRefs.back()).as<std::string>();
}
void property(const std::string &name, const std::string &value)
void property(const std::string_view name, const std::string &value)
{
mRow[name] = value;
for (auto &&[tag, ref] : mCachedRefs)
{
if (tag != name)
continue;
ref = value;
return;
}
mCachedRefs.emplace_back(name, mRow[name]);
std::get<1>(mCachedRefs.back()) = value;
}
int compare(const AtomImpl &b) const
......@@ -432,9 +439,11 @@ struct AtomImpl
Point mLocation;
int mRefcount;
cif::Row mRow;
mutable std::vector<std::tuple<std::string,cif::detail::ItemReference>> mCachedRefs;
mutable const Compound *mCompound = nullptr;
float mRadius = std::nanf("4");
mutable std::map<std::string, std::string> mCachedProperties;
bool mSymmetryCopy = false;
bool mClone = false;
......@@ -533,25 +542,25 @@ const cif::Row Atom::getRowAniso() const
}
template <>
std::string Atom::property<std::string>(const std::string &name) const
std::string Atom::property<std::string>(const std::string_view name) const
{
return impl()->property(name);
}
template <>
int Atom::property<int>(const std::string &name) const
int Atom::property<int>(const std::string_view name) const
{
auto v = impl()->property(name);
return v.empty() ? 0 : stoi(v);
}
template <>
float Atom::property<float>(const std::string &name) const
float Atom::property<float>(const std::string_view name) const
{
return stof(impl()->property(name));
}
void Atom::property(const std::string &name, const std::string &value)
void Atom::property(const std::string_view name, const std::string &value)
{
impl()->property(name, value);
}
......@@ -1736,7 +1745,7 @@ File::~File()
delete mImpl;
}
cif::Datablock& File::createDatablock(const std::string &name)
cif::Datablock& File::createDatablock(const std::string_view name)
{
auto db = new cif::Datablock(name);
......@@ -1978,6 +1987,58 @@ Atom Structure::getAtomByLabel(const std::string &atomID, const std::string &asy
throw std::out_of_range("Could not find atom with specified label");
}
Atom Structure::getAtomByPosition(Point p) const
{
double distance = std::numeric_limits<double>::max();
size_t index = std::numeric_limits<size_t>::max();
for (size_t i = 0; i < mAtoms.size(); ++i)
{
auto &a = mAtoms.at(i);
auto d = Distance(a.location(), p);
if (d < distance)
{
distance = d;
index = i;
}
}
if (index < mAtoms.size())
return mAtoms.at(index);
return {};
}
Atom Structure::getAtomByPositionAndType(Point p, std::string_view type, std::string_view res_type) const
{
double distance = std::numeric_limits<double>::max();
size_t index = std::numeric_limits<size_t>::max();
for (size_t i = 0; i < mAtoms.size(); ++i)
{
auto &a = mAtoms.at(i);
if (a.labelCompID() != res_type)
continue;
if (a.labelAtomID() != type)
continue;
auto d = Distance(a.location(), p);
if (d < distance)
{
distance = d;
index = i;
}
}
if (index < mAtoms.size())
return mAtoms.at(index);
return {};
}
const Residue &Structure::getResidue(const std::string &asymID, const std::string &compID, int seqID) const
{
for (auto &poly : mPolymers)
......@@ -2014,12 +2075,35 @@ const Residue &Structure::getResidue(const std::string &asymID, const std::strin
throw std::out_of_range("Could not find residue " + asymID + '/' + std::to_string(seqID));
}
Residue &Structure::getResidue(const std::string &asymID, const std::string &compID, int seqID)
{
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID, compID, seqID));
}
const Residue &Structure::getResidue(const std::string &asymID) const
{
for (auto &res : mNonPolymers)
{
if (res.asymID() != asymID)
continue;
return res;
}
throw std::out_of_range("Could not find residue " + asymID);
}
Residue &Structure::getResidue(const std::string &asymID)
{
return const_cast<Residue&>(const_cast<Structure const&>(*this).getResidue(asymID));
}
File &Structure::getFile() const
{
return mFile;
}
cif::Category &Structure::category(const char *name) const
cif::Category &Structure::category(std::string_view name) const
{
auto &db = datablock();
return db[name];
......@@ -2352,7 +2436,7 @@ void Structure::moveAtom(Atom &a, Point p)
a.location(p);
}
void Structure::changeResidue(const Residue &res, const std::string &newCompound,
void Structure::changeResidue(Residue &res, const std::string &newCompound,
const std::vector<std::tuple<std::string, std::string>> &remappedAtoms)
{
using namespace cif::literals;
......@@ -2416,6 +2500,8 @@ void Structure::changeResidue(const Residue &res, const std::string &newCompound
else
insertCompound(newCompound, false);
res.setCompoundID(newCompound);
auto &atomSites = db["atom_site"];
auto atoms = res.atoms();
......@@ -2437,8 +2523,15 @@ void Structure::changeResidue(const Residue &res, const std::string &newCompound
if (r.size() != 1)
continue;
if (a1 != a2)
r.front()["label_atom_id"] = a2;
if (a2.empty() or a2 == ".")
removeAtom(*i);
else if (a1 != a2)
{
auto ra = r.front();
ra["label_atom_id"] = a2;
ra["auth_atom_id"] = a2;
ra["type_symbol"] = AtomTypeTraits(compound->getAtomByID(a2).typeSymbol).symbol();
}
}
for (auto a : atoms)
......
src/Symmetry.cpp
......@@ -90,4 +90,66 @@ int GetSpacegroupNumber(std::string spacegroup)
return result;
}
// --------------------------------------------------------------------
int GetSpacegroupNumber(std::string spacegroup, SpacegroupName type)
{
if (spacegroup == "P 21 21 2 A")
spacegroup = "P 21 21 2 (a)";
else if (spacegroup.empty())
throw std::runtime_error("No spacegroup, cannot continue");
int result = 0;
if (type == SpacegroupName::full)
{
const size_t N = kNrOfSpaceGroups;
int32_t L = 0, R = static_cast<int32_t>(N - 1);
while (L <= R)
{
int32_t i = (L + R) / 2;
int d = spacegroup.compare(kSpaceGroups[i].name);
if (d > 0)
L = i + 1;
else if (d < 0)
R = i - 1;
else
{
result = kSpaceGroups[i].nr;
break;
}
}
}
else if (type == SpacegroupName::xHM)
{
for (auto &sg : kSpaceGroups)
{
if (sg.xHM == spacegroup)
{
result = sg.nr;
break;
}
}
}
else
{
for (auto &sg : kSpaceGroups)
{
if (sg.Hall == spacegroup)
{
result = sg.nr;
break;
}
}
}
// not found, see if we can find a match based on xHM name
if (result == 0)
throw std::runtime_error("Spacegroup name " + spacegroup + " was not found in table");
return result;
}
}
test/pdb2cif-test.cpp
......@@ -18,7 +18,6 @@ int main(int argc, char* argv[])
desc.add_options()
("input,i", po::value<std::string>(), "Input file")
("help,h", "Display help message")
("version", "Print version")
("verbose,v", "Verbose output")
("debug,d", po::value<int>(), "Debug level (for even more verbose output)");
......@@ -29,12 +28,6 @@ int main(int argc, char* argv[])
po::store(po::command_line_parser(argc, argv).options(desc).positional(p).run(), vm);
po::notify(vm);
if (vm.count("version"))
{
std::cout << argv[0] << " version " PACKAGE_VERSION << std::endl;
exit(0);
}
if (vm.count("help") or vm.count("input") == 0)
{
std::cerr << desc << std::endl;
......
test/unit-test.cpp
......@@ -33,6 +33,10 @@
// #include "cif++/DistanceMap.hpp"
#include "cif++/Cif++.hpp"
#include "cif++/BondMap.hpp"
#include "cif++/CifValidator.hpp"
namespace tt = boost::test_tools;
std::filesystem::path gTestDir = std::filesystem::current_path(); // filled in first test
......@@ -259,8 +263,10 @@ save__cat_2.desc
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -387,8 +393,10 @@ save__cat_1.c
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -535,8 +543,10 @@ save__cat_2.desc
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -741,8 +751,10 @@ save__cat_2.parent_id3
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -963,8 +975,10 @@ cat_2 3 cat_2:cat_1:3
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -1389,9 +1403,10 @@ cat_2 1 '_cat_2.num' '_cat_3.num' cat_3
} buffer(const_cast<char*>(dict), sizeof(dict) - 1);
std::istream is_dict(&buffer);
cif::Validator validator("test", is_dict);
cif::File f;
f.loadDictionary(is_dict);
f.setValidator(&validator);
// --------------------------------------------------------------------
......@@ -1687,4 +1702,58 @@ BOOST_AUTO_TEST_CASE(reading_file_1)
cif::File file;
BOOST_CHECK_THROW(file.load(is), std::runtime_error);
}
\ No newline at end of file
}
// 3d tests
using namespace mmcif;
BOOST_AUTO_TEST_CASE(t1)
{
// std::random_device rnd;
// std::mt19937 gen(rnd());
// std::uniform_real_distribution<float> dis(0, 1);
// Quaternion q{ dis(gen), dis(gen), dis(gen), dis(gen) };
// q = Normalize(q);
// Quaternion q{ 0.1, 0.2, 0.3, 0.4 };
Quaternion q{ 0.5, 0.5, 0.5, 0.5 };
q = Normalize(q);
const auto &&[angle0, axis0] = QuaternionToAngleAxis(q);
std::vector<Point> p1{
{ 16.979, 13.301, 44.555 },
{ 18.150, 13.525, 43.680 },
{ 18.656, 14.966, 43.784 },
{ 17.890, 15.889, 44.078 },
{ 17.678, 13.270, 42.255 },
{ 16.248, 13.734, 42.347 },
{ 15.762, 13.216, 43.724 }
};
auto p2 = p1;
Point c1 = CenterPoints(p1);
for (auto &p : p2)
p.rotate(q);
Point c2 = CenterPoints(p2);
auto q2 = AlignPoints(p1, p2);
const auto &&[angle, axis] = QuaternionToAngleAxis(q2);
BOOST_TEST(std::fmod(360 + angle, 360) == std::fmod(360 - angle0, 360), tt::tolerance(0.01));
for (auto &p : p1)
p.rotate(q2);
float rmsd = RMSd(p1, p2);
BOOST_TEST(rmsd < 1e-5);
// std::cout << "rmsd: " << RMSd(p1, p2) << std::endl;
}
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