source: src/Graph/BondGraph.cpp@ cd5aa0

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010-2012 University of Bonn. All rights reserved.
 * Copyright (C)  2013 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/>.
 */

/*
 * bondgraph.cpp
 *
 *  Created on: Oct 29, 2009
 *      Author: heber
 */

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

// boost::graph uses placement new
#include <boost/graph/adjacency_list.hpp>

#include "CodePatterns/MemDebug.hpp"

#include <algorithm>
#include <boost/bimap.hpp>
#include <boost/bind.hpp>
#include <boost/foreach.hpp>
#include <boost/function.hpp>
#include <boost/graph/max_cardinality_matching.hpp>
#include <deque>
#include <iostream>

#include "Atom/atom.hpp"
#include "Bond/bond.hpp"
#include "Graph/BondGraph.hpp"
#include "Box.hpp"
#include "Element/element.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Range.hpp"
#include "CodePatterns/Verbose.hpp"
#include "molecule.hpp"
#include "Element/periodentafel.hpp"
#include "Fragmentation/MatrixContainer.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "LinearAlgebra/Vector.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

const double BondGraph::BondThreshold = 0.4;   //!< CSD threshold in bond check which is the width of the interval whose center is the sum of the covalent radii

BondGraph::BondGraph() :
    BondLengthMatrix(NULL),
    IsAngstroem(true)
{}

BondGraph::BondGraph(bool IsA) :
    BondLengthMatrix(NULL),
    IsAngstroem(IsA)
{}

BondGraph::~BondGraph()
{
  CleanupBondLengthTable();
}

void BondGraph::CleanupBondLengthTable()
{
  if (BondLengthMatrix != NULL) {
    delete(BondLengthMatrix);
  }
}

bool BondGraph::LoadBondLengthTable(
    std::istream &input)
{
  Info FunctionInfo(__func__);
  bool status = true;
  MatrixContainer *TempContainer = NULL;

  // allocate MatrixContainer
  if (BondLengthMatrix != NULL) {
    LOG(1, "MatrixContainer for Bond length already present, removing.");
    delete(BondLengthMatrix);
    BondLengthMatrix = NULL;
  }
  TempContainer = new MatrixContainer;

  // parse in matrix
  if ((input.good()) && (status = TempContainer->ParseMatrix(input, 0, 1, 0))) {
    LOG(1, "Parsing bond length matrix successful.");
  } else {
    ELOG(1, "Parsing bond length matrix failed.");
    status = false;
  }

  if (status) // set to not NULL only if matrix was parsed
    BondLengthMatrix = TempContainer;
  else {
    BondLengthMatrix = NULL;
    delete(TempContainer);
  }
  return status;
}

double BondGraph::GetBondLength(
    int firstZ,
    int secondZ) const
{
  double return_length;
  if ((firstZ < 0) || (firstZ >= (int)BondLengthMatrix->Matrix[0].size()))
    return -1.;
  if ((secondZ < 0) || (secondZ >= (int)BondLengthMatrix->Matrix[0][firstZ].size()))
    return -1.;
  if (BondLengthMatrix == NULL) {
    return_length = -1.;
  } else {
    return_length = BondLengthMatrix->Matrix[0][firstZ][secondZ];
  }
  LOG(4, "INFO: Request for length between " << firstZ << " and "
      << secondZ << ": " << return_length << ".");
  return return_length;
}

range<double> BondGraph::CovalentMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = OtherWalker->getCovalentRadius() + Walker->getCovalentRadius();
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::BondLengthMatrixMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  ASSERT(BondLengthMatrix, "BondGraph::BondLengthMatrixMinMaxDistance() called without NULL BondLengthMatrix.");
  ASSERT(Walker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal element given.");
  ASSERT(OtherWalker, "BondGraph::BondLengthMatrixMinMaxDistance() - illegal other element given.");
  range<double> MinMaxDistance(0.,0.);
  MinMaxDistance.first = GetBondLength(Walker->getAtomicNumber()-1, OtherWalker->getAtomicNumber()-1);
  MinMaxDistance.first *= (IsAngstroem) ? 1. : 1. / AtomicLengthToAngstroem;
  MinMaxDistance.last = MinMaxDistance.first + BondThreshold;
  MinMaxDistance.first -= BondThreshold;
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const element * const Walker,
    const element * const OtherWalker) const
{
  range<double> MinMaxDistance(0.,0.);
  if (BondLengthMatrix == NULL) {// safety measure if no matrix has been parsed yet
    LOG(5, "DEBUG: Using Covalent radii criterion for [min,max) distances.");
    MinMaxDistance = CovalentMinMaxDistance(Walker, OtherWalker);
  } else {
    LOG(5, "DEBUG: Using tabulated bond table criterion for [min,max) distances.");
    MinMaxDistance = BondLengthMatrixMinMaxDistance(Walker, OtherWalker);
  }
  return MinMaxDistance;
}

range<double> BondGraph::getMinMaxDistance(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  return getMinMaxDistance(Walker->getType(), OtherWalker->getType());
}

range<double> BondGraph::getMinMaxDistanceSquared(
    const BondedParticle * const Walker,
    const BondedParticle * const OtherWalker) const
{
  // use non-squared version
  range<double> MinMaxDistance(getMinMaxDistance(Walker, OtherWalker));
  // and square
  MinMaxDistance.first *= MinMaxDistance.first;
  MinMaxDistance.last *= MinMaxDistance.last;
  return MinMaxDistance;
}

LinkedCell::LinkedCell_View BondGraph::getLinkedCell(const double max_distance) const
{
  return World::getInstance().getLinkedCell(max_distance);
}

std::set< const element *> BondGraph::getElementSetFromNumbers(const std::set<atomicNumber_t> &Set) const
{
  std::set< const element *> PresentElements;
  for(std::set<atomicNumber_t>::const_iterator iter = Set.begin(); iter != Set.end(); ++iter)
    PresentElements.insert( World::getInstance().getPeriode()->FindElement(*iter) );
  return PresentElements;
}

Box &BondGraph::getDomain() const
{
  return World::getInstance().getDomain();
}

unsigned int BondGraph::getTime() const
{
  return WorldTime::getTime();
}

bool BondGraph::operator==(const BondGraph &other) const
{
  if (IsAngstroem != other.IsAngstroem)
    return false;
  if (((BondLengthMatrix != NULL) && (other.BondLengthMatrix == NULL))
      || ((BondLengthMatrix == NULL) && (other.BondLengthMatrix != NULL)))
    return false;
  if ((BondLengthMatrix != NULL) && (other.BondLengthMatrix != NULL))
    if (*BondLengthMatrix != *other.BondLengthMatrix)
      return false;
  return true;
}

/** Corrects the bond degree by one at most if necessary.
 *
 * We are constraint to bonds in \a egdes, if the candidate bond is in edges, we
 * just don't increment its bond degree. We do not choose an alternative for this
 * atom.
 *
 * \param edges only these edges can be updated
 * \return number of corrections done
 */
int BondGraph::CorrectBondDegree(atom *_atom, const std::set<bond::ptr>& edges) const
{
  int NoBonds = 0;
  int OtherNoBonds = 0;
  int FalseBondDegree = 0;
  int CandidateDeltaNoBonds = -1;
  atom *OtherWalker = NULL;
  bond::ptr CandidateBond;

  NoBonds = _atom->CountBonds();
  LOG(3, "Walker " << *_atom << ": " << (int)_atom->getType()->getNoValenceOrbitals() << " > " << NoBonds << "?");
  const int DeltaNoBonds = (int)(_atom->getType()->getNoValenceOrbitals()) - NoBonds;
  if (DeltaNoBonds > 0) { // we have a mismatch, check all bonding partners for mismatch
    const BondList& ListOfBonds = _atom->getListOfBonds();
    for (BondList::const_iterator Runner = ListOfBonds.begin(); Runner != ListOfBonds.end(); (++Runner)) {
      OtherWalker = (*Runner)->GetOtherAtom(_atom);
      OtherNoBonds = OtherWalker->CountBonds();
      const int OtherDeltaNoBonds = (int)(OtherWalker->getType()->getNoValenceOrbitals()) - OtherNoBonds;
      LOG(3, "OtherWalker " << *OtherWalker << ": " << (int)OtherWalker->getType()->getNoValenceOrbitals() << " > " << OtherNoBonds << "?");
      if (OtherDeltaNoBonds > 0) { // check if possible candidate
        const BondList& OtherListOfBonds = OtherWalker->getListOfBonds();
        if ((CandidateBond == NULL)
            || (ListOfBonds.size() > OtherListOfBonds.size()) // pick the one with fewer number of bonds first
            || ((CandidateDeltaNoBonds < 0) || (CandidateDeltaNoBonds > OtherDeltaNoBonds)) ) // pick the one closer to fulfilling its valence first
        {
          CandidateDeltaNoBonds = OtherDeltaNoBonds;
          CandidateBond = (*Runner);
          LOG(3, "New candidate is " << *CandidateBond << ".");
        }
      }
    }
    if ((CandidateBond != NULL)) {
      if (edges.find(CandidateBond) != edges.end()) {
        CandidateBond->setDegree(CandidateBond->getDegree()+1);
        LOG(2, "Increased bond degree for bond " << *CandidateBond << ".");
      } else
        LOG(2, "Bond " << *CandidateBond << " is not in edges.");
    } else {
      ELOG(2, "Could not find correct degree for atom " << *_atom << ".");
      FalseBondDegree++;
    }
  }
  return FalseBondDegree;
}

/** Sets the weight of all connected bonds to one.
 */
void BondGraph::resetBondDegree(atom *_atom) const
{
  const BondList &ListOfBonds = _atom->getListOfBonds();
  for (BondList::const_iterator BondRunner = ListOfBonds.begin();
      BondRunner != ListOfBonds.end();
      ++BondRunner)
    (*BondRunner)->setDegree(1);
}

std::set<bond::ptr> BondGraph::getBackEdges() const
{
  DepthFirstSearchAnalysis DFS;
  DFS();

  const std::deque<bond::ptr>& backedgestack = DFS.getBackEdgeStack();
  std::set<bond::ptr> backedges(backedgestack.begin(), backedgestack.end());

  return backedges;
}

std::set<bond::ptr> BondGraph::getTreeEdges() const
{
  const std::set<bond::ptr> backedges = getBackEdges();
  std::set<bond::ptr> alledges;
  const World& world = World::getInstance();
  for(World::AtomConstIterator iter = world.getAtomIter();
      iter != world.atomEnd(); ++iter) {
    const BondList &ListOfBonds = (*iter)->getListOfBonds();
    alledges.insert(ListOfBonds.begin(), ListOfBonds.end());
  }
  std::set<bond::ptr> treeedges;
  std::set_difference(
      alledges.begin(), alledges.end(),
      backedges.begin(), backedges.end(),
      std::inserter(treeedges, treeedges.begin()));
  return treeedges;
}

struct hasDelta {
  bool operator()(atom *_atom) {
    const double delta =
        _atom->getType()->getNoValenceOrbitals() - _atom->CountBonds();
    return (fabs(delta) > 0);
  }
};

typedef std::set<atom *> deltaatoms_t;
typedef std::set<bond::ptr> deltaedges_t;

static int gatherAllDeltaAtoms(
    const deltaatoms_t &allatoms,
    deltaatoms_t &deltaatoms)
{
  static hasDelta delta;
  deltaatoms.clear();
  for (deltaatoms_t::const_iterator iter = allatoms.begin();
      iter != allatoms.end(); ++iter) {
    if (delta(*iter))
      deltaatoms.insert(*iter);
  }
  return deltaatoms.size();
}

typedef boost::bimap<int,atom*> AtomIndexLookup_t;

static AtomIndexLookup_t createAtomIndexLookup(
    const deltaedges_t &edges)
{
  AtomIndexLookup_t AtomIndexLookup;
  size_t index = 0;
  for (deltaedges_t::const_iterator edgeiter = edges.begin();
      edgeiter != edges.end(); ++edgeiter) {
    const bond::ptr & _bond = *edgeiter;

    // insert left into lookup
    std::pair< AtomIndexLookup_t::iterator, bool > inserter =
        AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->leftatom) );
    if (inserter.second)
      ++index;

    // insert right into lookup
    inserter = AtomIndexLookup.insert( AtomIndexLookup_t::value_type(index, _bond->rightatom) );
    if (inserter.second)
      ++index;
  }
  return AtomIndexLookup;
}

typedef std::map< std::pair<atom *, atom*>, bond::ptr> BondLookup_t;
static BondLookup_t createBondLookup(
    const deltaedges_t &edges)
{
  BondLookup_t BondLookup;
  for (deltaedges_t::const_iterator edgeiter = edges.begin();
      edgeiter != edges.end(); ++edgeiter) {
    const bond::ptr & _bond = *edgeiter;

    // insert bond into pair lookup
    BondLookup.insert(
        std::make_pair(
            std::make_pair(_bond->leftatom, _bond->rightatom),
            _bond)
    );
  }
  return BondLookup;
}

typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> graph_t;
typedef std::vector<boost::graph_traits<graph_t>::vertex_descriptor> mate_t;

static void fillEdgesIntoMatching(
    graph_t &g,
    mate_t &mate,
    const AtomIndexLookup_t &AtomIndexLookup,
    const BondLookup_t &BondLookup,
    deltaedges_t &matching
    )
{
  matching.clear();
  boost::graph_traits<graph_t>::vertex_iterator vi, vi_end;
  for(tie(vi,vi_end) = boost::vertices(g); vi != vi_end; ++vi)
    if (mate[*vi] != boost::graph_traits<graph_t>::null_vertex() && *vi < mate[*vi]) {
      atom * leftatom = AtomIndexLookup.left.at(*vi);
      atom * rightatom = AtomIndexLookup.left.at(mate[*vi]);
      std::pair<atom *,atom *> AtomPair(leftatom, rightatom);
      std::pair<atom *,atom *> reverseAtomPair(rightatom, leftatom);
      BondLookup_t::const_iterator iter = BondLookup.find(AtomPair);
      if (iter != BondLookup.end()) {
        matching.insert(iter->second);
      } else {
        iter = BondLookup.find(reverseAtomPair);
        if (iter != BondLookup.end()) {
          matching.insert(iter->second);
        } else
          ASSERT(0, "fillEdgesIntoMatching() - cannot find ("+toString(*vi)+","
              +toString(mate[*vi])+") in BondLookup.");
      }
    }
}

static bool createMatching(
    deltaedges_t &deltaedges,
    deltaedges_t &matching
    )
{
  // gather all vertices for graph in a lookup structure
  // to translate the graphs indices to atoms and also to associated bonds
  AtomIndexLookup_t AtomIndexLookup = createAtomIndexLookup(deltaedges);
  BondLookup_t BondLookup = createBondLookup(deltaedges);
  const int n_vertices = AtomIndexLookup.size();

  // construct graph
  graph_t g(n_vertices);

  // add edges to graph
  for (deltaedges_t::const_iterator edgeiter = deltaedges.begin();
      edgeiter != deltaedges.end(); ++edgeiter) {
    const bond::ptr & _bond = *edgeiter;
    boost::add_edge(
        AtomIndexLookup.right.at(_bond->leftatom),
        AtomIndexLookup.right.at(_bond->rightatom),
        g);
  }

  // mate structure contains unique edge partner to every vertex in matching
  mate_t mate(n_vertices);

  // get maximum matching over given edges
  bool success = boost::checked_edmonds_maximum_cardinality_matching(g, &mate[0]);

  if (success) {
    LOG(1, "STATUS: Found a matching of size " << matching_size(g, &mate[0]) << ".");
    fillEdgesIntoMatching(g, mate, AtomIndexLookup, BondLookup, matching);
  } else {
    ELOG(2, "Failed to find maximum matching.");
  }

  return success;
}


int BondGraph::calculateBondDegree(const deltaatoms_t &allatoms) const
{
  deltaatoms_t deltaatoms;
  deltaedges_t deltaedges;
  deltaedges_t matching;

  LOG(1, "STATUS: Calculating bond degrees.");
  if (DoLog(2)) {
    std::stringstream output;
    output << "INFO: All atoms are: ";
    BOOST_FOREACH( atom *_atom, allatoms) {
      output << *_atom << " ";
    }
    LOG(2, output.str());
  }

  size_t IterCount = 0;
  // 1. gather all atoms with delta > 0
  while ((gatherAllDeltaAtoms(allatoms, deltaatoms) != 0) && (IterCount < 100)) {
    if (DoLog(3)) {
      std::stringstream output;
      output << "DEBUG: Current delta atoms are: ";
      BOOST_FOREACH( atom *_atom, deltaatoms) {
        output << *_atom << " ";
      }
      LOG(3, output.str());
    }

    // 2. gather all edges that connect atoms with both(!) having delta > 0
    deltaedges.clear();
    for (deltaatoms_t::const_iterator atomiter = deltaatoms.begin();
        atomiter != deltaatoms.end(); ++atomiter) {
      const atom * const _atom = *atomiter;
      const BondList& ListOfBonds = (*atomiter)->getListOfBonds();
      for (BondList::const_iterator bonditer = ListOfBonds.begin();
          bonditer != ListOfBonds.end();++bonditer) {
        const bond::ptr _bond = *bonditer;
        if ((_bond->leftatom == _atom) && (deltaatoms.find(_bond->rightatom) != deltaatoms.end()))
          deltaedges.insert(_bond);
        else if ((_bond->rightatom == _atom) && (deltaatoms.find(_bond->leftatom) != deltaatoms.end()))
          deltaedges.insert(_bond);
      }
    }
    if (DoLog(3)) {
      std::stringstream output;
      output << "DEBUG: Current delta bonds are: ";
      BOOST_FOREACH( bond::ptr _bond, deltaedges) {
        output << *_bond << " ";
      }
      LOG(3, output.str());
    }
    if (deltaedges.empty())
      break;

    // 3. build matching over these edges
    if (!createMatching(deltaedges, matching))
      break;
    if (DoLog(2)) {
      std::stringstream output;
      output << "DEBUG: Resulting matching is: ";
      BOOST_FOREACH( bond::ptr _bond, matching) {
        output << *_bond << " ";
      }
      LOG(2, output.str());
    }

    // 4. increase degree for each edge of the matching
    LOG(2, "DEBUG: Increasing bond degrees by one.");
    for (deltaedges_t::iterator edgeiter = matching.begin();
        edgeiter != matching.end(); ++edgeiter) {
      (*edgeiter)->setDegree( (*edgeiter)->getDegree()+1 );
    }

    // 6. stop after a given number of iterations or when all atoms are happy.
    ++IterCount;
  };

  // check whether we still have deltaatoms
  {
    hasDelta delta;
    for (deltaatoms_t::const_iterator iter = allatoms.begin();
        iter != allatoms.end(); ++iter)
      if (delta(*iter))
        ELOG(2, "Could not find satisfy charge neutrality for atom " << **iter << ".");
  }

  return deltaedges.size();
}