source: src/Graph/Graph6Writer.cpp@ 1c0b0b

/*
 * 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 <http://www.gnu.org/licenses/>.
 */

/*
 * Graph6Writer.cpp
 *
 *  Created on: Apr 2, 2021
 *      Author: heber
 */


// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "Graph6Writer.hpp"

#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"

#include <cassert>
#include <cmath>
#include <iostream>

#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<const atom *> 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<<n_pos))>>((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; j<n; ++j)
    for (size_t i=0; i<j; ++i) {
      // std::cout << "\t\n" << (int)value << " " << byte_pos << std::endl;

      unsigned int bit = _atoms[i]->IsBondedTo(_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 atom *>::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<atomId_t> gathered_atoms = gatherer(start_atom_id);

  // print all nodes
  const World& world = World::getConstInstance();
  for (std::vector<atomId_t>::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();
  }
}