/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2021 Frederik Heber. All rights reserved.
*
*
* This file is part of MoleCuilder.
*
* MoleCuilder is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* MoleCuilder is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MoleCuilder. If not, see .
*/
/*
* Graph6Writer.cpp
*
* Created on: Apr 2, 2021
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "Graph6Writer.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include
#include
#include
#include "Atom/atom.hpp"
#include "Descriptors/AtomIdDescriptor.hpp"
#include "Element/element.hpp"
#include "Graph/BoostGraphCreator.hpp"
#include "Graph/BreadthFirstSearchGatherer.hpp"
#include "World.hpp"
//#include "CodePatterns/MemDebug.hpp"
Graph6Writer::Graph6Writer(const std::vector atoms):
_atoms(atoms)
{}
void Graph6Writer::write_n(std::ostream& out) {
const unsigned long n = _atoms.size();
if (n<62) {
out << ((unsigned char)(n+63));
return;
}
out << ((unsigned char)126);
int num_bytes = 2;
if (n> 258047) {
out << ((unsigned char)126);
num_bytes = 3;
}
for(int value=num_bytes; value>=0; value--) {
unsigned char c = 0;
int n_pos = 6*(value+1)-1;
for(int c_pos=5; c_pos>=0; n_pos--, c_pos--) {
c += (n & (1<>((int)n_pos/6);
}
out << (c+63);
}
}
/* Given an iterator over the adjacency matrix in the order (0,1),(0,2),(1,2),(0,3),(1,3),(2,3),...,(n-1,n)
this writes a graph6 representation to out. */
void Graph6Writer::write_graph6(std::ostream& out) {
write_n(out);
const unsigned long n = _atoms.size();
unsigned char value = 0;
int byte_pos = 5;
unsigned int bytes_written = 0;
for (size_t j=0; jIsBondedTo(_atoms[j]);
LOG(2, "DEBUG: (" << i << "," << j << ") = " << bit << "," << value << " | " << bit << " << " << byte_pos << " = " << (unsigned int)value << " | " << (bit << byte_pos));
value = value | (bit << byte_pos--);
if (byte_pos < 0) {
LOG(2, "DEBUG: Writing byte " << value << " into range [" << (unsigned char)63 << "," << (unsigned char)126 << "]");
ASSERT( (value+63) <= 126,
"Graph6Writer::write_graph6() - char to write is outside "+toString((unsigned char)63)
+" and "+toString((unsigned char)126));
out << (unsigned char)(value+63);
bytes_written++;
value = 0;
byte_pos = 5;
}
}
if (byte_pos!=5) {
ASSERT( (value+63) <= 126,
"Graph6Writer::write_graph6() - char to write is outside "+toString((unsigned char)63)
+" and "+toString((unsigned char)126));
LOG(2, "DEBUG: Writing byte " << value << " into range [" << (unsigned char)63 << "," << (unsigned char)126 << "]");
out << (unsigned char)(value+63);
bytes_written++;
value=0;
}
ASSERT( value==0,
"Graph6Writer::write_graph6() - byte is not null, i.e. chars left to write?");
ASSERT( bytes_written == (unsigned int)ceil(n*(n-1)/12.0f),
"Graph6Writer::write_graph6() - unexpected number of bytes written");
}
/**
* Picks a non-hydrogen from all hydrogen atoms in the current set
* of atoms.
*
* Returns -1 if none could be found.
*/
atomId_t Graph6Writer::getBoundaryNonHydrogen() const {
atomId_t start_atom_id = -1;
int bond_degree = 16;
for(std::vector::const_iterator iter = _atoms.begin();
iter != _atoms.end(); ++iter) {
const atom *walker = *iter;
if (walker->getElement().getSymbol() == "H") {
const BondList& bond_list = walker->getListOfBonds();
assert(bond_list.size() == 1);
const atom *other = bond_list.front().get()->GetOtherAtom(walker);
const int number_of_bonds = other->getListOfBonds().size();
if ((other->getElement().getSymbol() != "H") && (bond_degree > number_of_bonds)) {
start_atom_id = other->getId();
bond_degree = number_of_bonds;
}
}
}
return start_atom_id;
}
bool OnlyNonHydrogens(const bond &_bond) {
return _bond.HydrogenBond == 0;
}
void Graph6Writer::write_elementlist(std::ostream& out) {
/** Execute a Breadth-First Search discovery from one terminal atom (e.g.,
* pick random hydrogen and then it's bond-neighbor if it is non-hydrogen).
* Then return the element list in that ordering.
*
* The graph6 string does not account for the inherent graph symmetries
* (e.g., BW having 123<->321 but not 123<->132 symmetry).
*/
// pick bond neighbor of a hydrogen atom
atomId_t start_atom_id = getBoundaryNonHydrogen();
if (start_atom_id == (unsigned int)-1) {
// fall back to first atom in list
start_atom_id = _atoms.front()->getId();
}
// do an unlimited BFS and get set of nodes, ordered by discovery level
BoostGraphCreator graphCreator;
graphCreator.createFromAtoms(_atoms, OnlyNonHydrogens);
BreadthFirstSearchGatherer gatherer(graphCreator);
const std::vector gathered_atoms = gatherer(start_atom_id);
// print all nodes
const World& world = World::getConstInstance();
for (std::vector::const_iterator iter = gathered_atoms.begin();
iter != gathered_atoms.end(); ++iter) {
if (iter != gathered_atoms.begin())
out << ' ';
const atom* walker = world.getAtom(AtomById(*iter));
assert(walker != NULL);
out << walker->getElement().getSymbol();
}
}