Skip to content

L
libcifpp
Commit 9b1e9356 by maarten
Options

met zonder resources, optioneel 

git-svn-id: svn+ssh://gitlab/srv/svn-repos/pdb-redo/trunk@475 a1961a4f-ab94-4bcc-80e8-33b5a54de466
parent f1dfe12c
Showing with 232 additions and 107 deletions
include/cif++/Cif++.h
......@@ -15,8 +15,6 @@
#include "cif++/CifUtils.h"
extern int VERBOSE;
/*
Simple C++ interface to CIF files.
......@@ -83,6 +81,9 @@ extern int VERBOSE;
namespace cif
{
// flag for verbose output
extern int VERBOSE;
using std::string;
using std::vector;
......
include/cif++/mrsrc.h
......@@ -10,6 +10,27 @@
#include <list>
#include <exception>
namespace mrsrc
{
class rsrc_not_found_exception : public std::exception
{
public:
virtual const char* what() const throw() { return "resource not found"; }
};
class rsrc;
typedef std::list<rsrc> rsrc_list;
}
// --------------------------------------------------------------------
// By default, we assume mrc is used to create the resources. If you need
// local disk storage instead, use the NO_RSRC macro. See below for
// usage.
#if defined(USE_RSRC)
/*
Resources are data sources for the application.
......@@ -47,15 +68,6 @@ extern const char gResourceName[];
namespace mrsrc
{
class rsrc_not_found_exception : public std::exception
{
public:
virtual const char* what() const throw() { return "resource not found"; }
};
class rsrc;
typedef std::list<rsrc> rsrc_list;
class rsrc
{
public:
......@@ -149,4 +161,121 @@ rsrc::rsrc(const std::string& path)
}
#else
// --------------------------------------------------------------------
// Fall back option for resources, locate the data in the path specified
// in then environment variable RESOURCE_DIR
#include <cstdlib>
#include <unistd.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <map>
namespace mrsrc
{
namespace detail
{
class rsrc_loader
{
public:
static rsrc_loader& instance()
{
static rsrc_loader sInstance;
return sInstance;
}
const std::vector<char>& load(const std::string& path)
{
if (not m_loaded.count(path))
{
std::string p = m_rsrc_dir + '/' + path;
std::ifstream file(p);
std::vector<char> result;
if (file.is_open())
{
std::streambuf* b = file.rdbuf();
size_t length = b->pubseekoff(0, std::ios::end);
b->pubseekpos(0);
result.resize(length);
b->sgetn(result.data(), length);
}
m_loaded.emplace(path, std::move(result));
}
return m_loaded.at(path);
}
private:
rsrc_loader()
{
const char* rsrc_dir = getenv("RESOURCE_DIR");
#if defined(RSRC_DIR)
if (rsrc_dir == nullptr)
rsrc_dir = RSRC_DIR;
#endif
if (rsrc_dir == nullptr)
{
char pb[PATH_MAX];
const char* cwd = getcwd(pb, PATH_MAX);
if (cwd == nullptr)
throw std::runtime_error("Could not locate resource directory nor current directory");
std::cerr << "RESOURCE_DIR not defined, falling back to " << cwd << std::endl;
rsrc_dir = cwd;
}
m_rsrc_dir = rsrc_dir;
}
std::string m_rsrc_dir;
std::map<std::string,std::vector<char>> m_loaded;
};
}
class rsrc
{
public:
rsrc(const rsrc&) = delete;
rsrc& operator=(const rsrc&) = delete;
rsrc(const std::string& path)
: m_data(detail::rsrc_loader::instance().load(path))
// , m_name(path)
{}
// std::string name() const { return m_name; }
const char* data() const { return m_data.data(); }
unsigned long size() const { return m_data.size(); }
explicit operator bool () const { return not m_data.empty(); }
// rsrc_list children() const;
private:
const std::vector<char>& m_data;
// const std::string m_name;
};
}
#endif
#endif
src/BondMap.cpp
......@@ -57,7 +57,7 @@ BondMap::BondMap(const Structure& p)
if (compounds.count(c))
continue;
if (VERBOSE > 1)
if (cif::VERBOSE > 1)
cerr << "Warning: mon_id " << c << " is missing in the chem_comp category" << endl;
compounds.insert(c);
}
......@@ -151,14 +151,14 @@ BondMap::BondMap(const Structure& p)
auto* compound = mmcif::Compound::create(c);
if (not compound)
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Missing compound information for " << c << endl;
continue;
}
if (compound->isWater())
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "skipping water in bond map calculation" << endl;
continue;
}
......
src/Cif++.cpp
......@@ -16,9 +16,7 @@
#include <boost/iostreams/filter/bzip2.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#if defined(USE_RSRC)
#include "cif++/mrsrc.h"
#endif
#include "cif++/Cif++.h"
#include "cif++/CifParser.h"
......@@ -30,11 +28,11 @@ namespace ba = boost::algorithm;
namespace fs = boost::filesystem;
namespace io = boost::iostreams;
extern int VERBOSE;
namespace cif
{
int VERBOSE = 0;
static const char* kEmptyResult = "";
// --------------------------------------------------------------------
......@@ -2690,7 +2688,6 @@ void File::loadDictionary(const char* dict)
}
catch (...) {}
#if defined(USE_RSRC)
mrsrc::rsrc dictData(dictFile.string());
if (dictData)
......@@ -2708,7 +2705,6 @@ void File::loadDictionary(const char* dict)
loadDictionary(is);
break;
}
#endif
throw runtime_error("Dictionary not found or defined (" + name + ")");
}
......
src/Cif2PDB.cpp
......@@ -2932,7 +2932,7 @@ int WriteHeterogen(ostream& pdbFile, Datablock& db)
h->numHetAtoms += 1;
}
if (VERBOSE > 1 and not missingHetNames.empty())
if (cif::VERBOSE > 1 and not missingHetNames.empty())
cerr << "Missing het name(s) for " << ba::join(missingHetNames, ", ") << endl;
boost::format kHET("HET %3.3s %1.1s%4.4d%1.1s %5.5d");
......
src/Compound.cpp
......@@ -15,9 +15,7 @@
#include "cif++/Compound.h"
#include "cif++/CifUtils.h"
#if defined(USE_RSRC)
#include "cif++/mrsrc.h"
#endif
using namespace std;
namespace ba = boost::algorithm;
......@@ -84,7 +82,7 @@ class IsomerDB
IsomerDB::IsomerDB()
{
#if defined(USE_RSRC)
// #if defined(USE_RSRC)
mrsrc::rsrc isomers("isomers.txt");
// mrsrc::rsrc isomers("isomers-with-sugar.xml");
if (not isomers)
......@@ -96,16 +94,16 @@ IsomerDB::IsomerDB()
} buffer(const_cast<char*>(isomers.data()), isomers.size());
istream is(&buffer);
#else
cerr << "resource support was not compiled in, falling back to a local file" << endl;
fs::path isomersFile = "isomers.txt";
if (not fs::exists(isomersFile))
isomersFile = fs::path(cif::get_executable_path()).parent_path() / "isomers.txt";
// #else
// cerr << "resource support was not compiled in, falling back to a local file" << endl;
// fs::path isomersFile = "isomers.txt";
// if (not fs::exists(isomersFile))
// isomersFile = fs::path(cif::get_executable_path()).parent_path() / "isomers.txt";
fs::ifstream is(isomersFile);
if (not is.is_open())
throw runtime_error("Could not open the file isomers.txt");
#endif
// fs::ifstream is(isomersFile);
// if (not is.is_open())
// throw runtime_error("Could not open the file isomers.txt");
// #endif
string line;
......@@ -357,7 +355,7 @@ Compound::Compound(const fs::path& file, const std::string& id,
b.type = delocalizedBond;
else
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Unimplemented chem_comp_bond.type " << type << " in " << id << endl;
b.type = singleBond;
}
......@@ -407,7 +405,7 @@ Compound::Compound(const fs::path& file, const std::string& id,
cc.volumeSign = both;
else
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Unimplemented chem_comp_chir.volume_sign " << volumeSign << " in " << id << endl;
continue;
}
......@@ -667,7 +665,7 @@ bool Compound::isIsomerOf(const Compound& c) const
vector<tuple<string,string>> mapping;
result = StructuresAreIsomeric(mAtoms, mBonds, c.mAtoms, c.mBonds, mapping);
if (VERBOSE and result)
if (cif::VERBOSE and result)
{
for (auto& m: mapping)
cerr << " " << get<0>(m) << " => " << get<1>(m) << endl;
......@@ -801,7 +799,7 @@ Link::Link(cif::Datablock& db)
b.type = delocalizedBond;
else
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Unimplemented chem_link_bond.type " << type << " in " << mId << endl;
b.type = singleBond;
}
......@@ -863,7 +861,7 @@ Link::Link(cif::Datablock& db)
cc.volumeSign = both;
else
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Unimplemented chem_link_chir.volume_sign " << volumeSign << " in " << mId << endl;
continue;
}
......
src/DistanceMap.cpp
......@@ -115,7 +115,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
sort(c.begin(), c.end());
float mz = median();
if (VERBOSE > 1)
if (cif::VERBOSE > 1)
cerr << "median position of atoms: " << Point(mx, my, mz) << endl;
auto calculateD = [&](float m, float c)
......@@ -136,7 +136,7 @@ DistanceMap::DistanceMap(const Structure& p, const clipper::Spacegroup& spacegro
if (mD.mX != 0 or mD.mY != 0 or mD.mZ != 0)
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "moving coorinates by " << mD.mX << ", " << mD.mY << " and " << mD.mZ << endl;
for_each(locations.begin(), locations.end(), [&](auto& p) { p += mD; });
......
src/PDB2CifRemark3.cpp
#include "cif++/Config.h"
#include "cif++/Cif++.h"
#include <map>
#include <set>
......@@ -996,7 +997,7 @@ string Remark3Parser::nextLine()
break;
}
if (VERBOSE >= 2)
if (cif::VERBOSE >= 2)
cerr << "RM3: " << mLine << endl;
return mLine;
......@@ -1010,7 +1011,7 @@ bool Remark3Parser::match(const char* expr, int nextState)
if (result)
mState = nextState;
else if (VERBOSE >= 3)
else if (cif::VERBOSE >= 3)
cerr << cif::coloured("No match:", cif::scWHITE, cif::scRED) << " '" << expr << '\'' << endl;
return result;
......@@ -1070,7 +1071,7 @@ float Remark3Parser::parse()
continue;
}
if (VERBOSE >= 2)
if (cif::VERBOSE >= 2)
cerr << cif::coloured("Dropping line:", cif::scWHITE, cif::scRED) << " '" << mLine << '\'' << endl;
++dropped;
......@@ -1122,7 +1123,7 @@ void Remark3Parser::storeCapture(const char* category, initializer_list<const ch
if (iequals(value, "NULL") or iequals(value, "NONE") or iequals(value, "Inf") or iequals(value, "+Inf") or iequals(value, string(value.length(), '*')))
continue;
if (VERBOSE >= 3)
if (cif::VERBOSE >= 3)
cerr << "storing: '" << value << "' in _" << category << '.' << item << endl;
auto& cat = mDb[category];
......@@ -1291,7 +1292,7 @@ bool Remark3Parser::parse(const string& expMethod, PDBRecord* r, cif::Datablock&
if (line != "REFINEMENT.")
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Unexpected data in REMARK 3" << endl;
return false;
}
......@@ -1303,7 +1304,7 @@ bool Remark3Parser::parse(const string& expMethod, PDBRecord* r, cif::Datablock&
if (not regex_match(line, m, rxp))
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Expected valid PROGRAM line in REMARK 3" << endl;
return false;
}
......@@ -1343,7 +1344,7 @@ bool Remark3Parser::parse(const string& expMethod, PDBRecord* r, cif::Datablock&
score = 0;
}
if (VERBOSE >= 2)
if (cif::VERBOSE >= 2)
cerr << "Score for " << parser->program() << ": " << score << endl;
if (score > 0)
......@@ -1382,7 +1383,7 @@ bool Remark3Parser::parse(const string& expMethod, PDBRecord* r, cif::Datablock&
tryParser(new TNT_Remark3Parser(program, expMethod, r, db));
else if (ba::starts_with(program, "X-PLOR"))
tryParser(new XPLOR_Remark3Parser(program, expMethod, r, db));
else if (VERBOSE)
else if (cif::VERBOSE)
cerr << "Skipping unknown program (" << program << ") in REMARK 3" << endl;
}
......@@ -1415,7 +1416,7 @@ bool Remark3Parser::parse(const string& expMethod, PDBRecord* r, cif::Datablock&
auto& best = scores.front();
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Choosing " << best.parser->program() << " version '" << best.parser->version() << "' as refinement program. Score = " << best.score << endl;
auto& software = db["software"];
......
src/Secondary.cpp
......@@ -753,7 +753,7 @@ void CalculateSecondaryStructure(Structure& s)
CalculateBetaSheets(residues);
CalculateAlphaHelices(residues);
if (VERBOSE)
if (cif::VERBOSE)
{
for (auto& r: residues)
{
......@@ -813,7 +813,7 @@ DSSPImpl::DSSPImpl(const Structure& s)
: mStructure(s)
, mPolymers(mStructure.polymers())
{
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Calculating DSSP ";
size_t nRes = accumulate(mPolymers.begin(), mPolymers.end(),
......@@ -835,18 +835,18 @@ DSSPImpl::DSSPImpl(const Structure& s)
mResidues[i + 1].assignHydrogen();
}
if (VERBOSE) cerr << ".";
if (cif::VERBOSE) cerr << ".";
CalculateHBondEnergies(mResidues);
if (VERBOSE) cerr << ".";
if (cif::VERBOSE) cerr << ".";
CalculateBetaSheets(mResidues);
if (VERBOSE) cerr << ".";
if (cif::VERBOSE) cerr << ".";
CalculateAlphaHelices(mResidues);
if (VERBOSE) cerr << endl;
if (cif::VERBOSE) cerr << endl;
if (VERBOSE > 1)
if (cif::VERBOSE > 1)
{
for (auto& r: mResidues)
{
......@@ -892,7 +892,7 @@ SecondaryStructureType DSSP::operator()(const string& inAsymID, int inSeqID) con
[&](auto& r) { return r.mM.asymID() == inAsymID and r.mM.seqID() == inSeqID; });
if (i != mImpl->mResidues.end())
result = i->mSecondaryStructure;
else if (VERBOSE)
else if (cif::VERBOSE)
cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << endl;
return result;
}
......@@ -916,7 +916,7 @@ bool DSSP::isAlphaHelixEndBeforeStart(const string& inAsymID, int inSeqID) const
if (i != mImpl->mResidues.end() and i + 1 != mImpl->mResidues.end())
result = i->GetHelixFlag(4) == helixEnd and (i + 1)->GetHelixFlag(4) == helixStart;
else if (VERBOSE)
else if (cif::VERBOSE)
cerr << "Could not find secondary structure for " << inAsymID << ':' << inSeqID << endl;
return result;
......
src/Statistics.cpp
......@@ -485,7 +485,7 @@ void StatsCollector::initialize()
float radius = shape.radius();
if (VERBOSE > 2)
if (cif::VERBOSE > 2)
cerr << (atomData.size() + 1) << '\t'
<< AtomTypeTraits(atom.type()).symbol() << '\t'
<< radius << endl;
......@@ -551,7 +551,7 @@ void StatsCollector::initialize()
double qa = dd * (swxs * swy - swx * swxy);
double qb = dd * (sw * swxy - swx * swy);
if (VERBOSE > 1)
if (cif::VERBOSE > 1)
{
swys = dd * (swys - (qa * swy + qb * swxy)) / (ns - 2);
cerr << endl
......@@ -560,7 +560,7 @@ void StatsCollector::initialize()
qb += 1.0;
if (VERBOSE > 1)
if (cif::VERBOSE > 1)
{
cerr << endl
<< "Rescale SD(delta-rho) using Q-Q plot for asym " << zsc.first << ':' << endl
......@@ -656,7 +656,7 @@ vector<ResidueStatistics> StatsCollector::collect(const vector<tuple<string,int,
float radius = shape.radius();
if (VERBOSE > 2)
if (cif::VERBOSE > 2)
cerr << (atomData.size() + 1) << '\t'
<< AtomTypeTraits(atom.type()).symbol() << '\t'
<< radius << endl;
......@@ -716,7 +716,7 @@ vector<ResidueStatistics> StatsCollector::collect(const vector<tuple<string,int,
{
if (compAtom.id == "OXT")
--n;
else if (VERBOSE > 1)
else if (cif::VERBOSE > 1)
cerr << "Missing atom '" << compAtom.id << "' in residue " << asymID << ':' << seqID << endl;
continue;
}
......@@ -794,7 +794,7 @@ ResidueStatistics StatsCollector::collect(const vector<Atom>& atoms) const
float radius = shape.radius();
if (VERBOSE > 2)
if (cif::VERBOSE > 2)
cerr << (atomData.size() + 1) << '\t'
<< AtomTypeTraits(atom.type()).symbol() << '\t'
<< radius << endl;
......@@ -995,7 +995,7 @@ EDIAStatsCollector::EDIAStatsCollector(MapMaker<float>& mm,
else
ediaBFactor = kAverageBFactors[i];
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Calculating radii with B Factor " << ediaBFactor << endl;
for (auto atom: mStructure.atoms())
......@@ -1006,7 +1006,7 @@ EDIAStatsCollector::EDIAStatsCollector(MapMaker<float>& mm,
AtomShape shape(atom, mResHigh, mResLow, mElectronScattering, ediaBFactor);
mRadii[atom.type()] = shape.radius();
if (VERBOSE)
if (cif::VERBOSE)
cerr << "Radius for atom with type " << AtomTypeTraits(atom.type()).symbol() << " is " << mRadii[atom.type()] << endl;
}
}
......@@ -1034,7 +1034,7 @@ void EDIAStatsCollector::calculate(vector<AtomData>& atomData) const
auto& atom = data.atom;
float radius = mRadii.at(atom.type());
// if (VERBOSE > 2)
// if (cif::VERBOSE > 2)
// cerr << (atomData.size() + 1) << '\t'
// << AtomTypeTraits(atom.type()).symbol() << '\t'
// << radius << endl;
......@@ -1127,7 +1127,7 @@ void EDIAStatsCollector::calculate(vector<AtomData>& atomData) const
}
}
// if (VERBOSE > 2)
// if (cif::VERBOSE > 2)
// cout << Point(p) << ":\td: " << xmap[iw] << "\tz: " << z << "\to: " << o << "\tzraw: " << ((xmap[iw] - meanDensity) / rmsDensity) << "\twp: " << wp << endl;
ediaSum[0] += z * wp * o;
......
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