source: src/molecule_graph.cpp@ d09093

/*
 * Project: MoleCuilder
 * Description: creates and alters molecular systems
 * Copyright (C)  2010 University of Bonn. All rights reserved.
 * Please see the LICENSE file or "Copyright notice" in builder.cpp for details.
 */

/*
 * molecule_graph.cpp
 *
 *  Created on: Oct 5, 2009
 *      Author: heber
 */

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

#include "CodePatterns/MemDebug.hpp"

#include <stack>

#include "atom.hpp"
#include "Bond/bond.hpp"
#include "Box.hpp"
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Info.hpp"
#include "CodePatterns/Log.hpp"
#include "CodePatterns/Verbose.hpp"
#include "config.hpp"
#include "Graph/DepthFirstSearchAnalysis.hpp"
#include "element.hpp"
#include "Graph/BondGraph.hpp"
#include "Helpers/defs.hpp"
#include "Helpers/fast_functions.hpp"
#include "Helpers/helpers.hpp"
#include "LinearAlgebra/RealSpaceMatrix.hpp"
#include "linkedcell.hpp"
#include "molecule.hpp"
#include "PointCloudAdaptor.hpp"
#include "World.hpp"
#include "WorldTime.hpp"

#define MAXBONDS 8


/** Fills the bond structure of this chain list subgraphs that are derived from a complete \a *reference molecule.
 * Calls this routine in each MoleculeLeafClass::next subgraph if it's not NULL.
 * \param *reference reference molecule with the bond structure to be copied
 * \param **&ListOfLocalAtoms Lookup table for this subgraph and index of each atom in \a *reference, may be NULL on start, then it is filled
 * \param FreeList true - ***ListOfLocalAtoms is free'd before return, false - it is not
 * \return true - success, false - failure
 */
bool molecule::FillBondStructureFromReference(const molecule * const reference, atom **&ListOfLocalAtoms, bool FreeList)
{
  atom *OtherWalker = NULL;
  atom *Father = NULL;
  bool status = true;
  int AtomNo;

  DoLog(1) && (Log() << Verbose(1) << "Begin of FillBondStructureFromReference." << endl);
  // fill ListOfLocalAtoms if NULL was given
  if (!FillListOfLocalAtoms(ListOfLocalAtoms, reference->getAtomCount())) {
    DoLog(1) && (Log() << Verbose(1) << "Filling of ListOfLocalAtoms failed." << endl);
    return false;
  }

  if (status) {
    DoLog(1) && (Log() << Verbose(1) << "Creating adjacency list for molecule " << getName() << "." << endl);
    // remove every bond from the list
    for_each(begin(), end(),
        boost::bind(&BondedParticle::ClearBondsAtStep,_1,WorldTime::getTime()));


    for(molecule::const_iterator iter = begin(); iter != end(); ++iter) {
      Father = (*iter)->GetTrueFather();
      AtomNo = Father->getNr(); // global id of the current walker
      const BondList& ListOfBonds = Father->getListOfBonds();
      for (BondList::const_iterator Runner = ListOfBonds.begin();
          Runner != ListOfBonds.end();
          ++Runner) {
        OtherWalker = ListOfLocalAtoms[(*Runner)->GetOtherAtom((*iter)->GetTrueFather())->getNr()]; // local copy of current bond partner of walker
        if (OtherWalker != NULL) {
          if (OtherWalker->getNr() > (*iter)->getNr())
            AddBond((*iter), OtherWalker, (*Runner)->BondDegree);
        } else {
          DoLog(1) && (Log() << Verbose(1) << "OtherWalker = ListOfLocalAtoms[" << (*Runner)->GetOtherAtom((*iter)->GetTrueFather())->getNr() << "] is NULL!" << endl);
          status = false;
        }
      }
    }
  }

  if ((FreeList) && (ListOfLocalAtoms != NULL)) {
    // free the index lookup list
    delete[](ListOfLocalAtoms);
  }
  DoLog(1) && (Log() << Verbose(1) << "End of FillBondStructureFromReference." << endl);
  return status;
};

/** Checks for presence of bonds within atom list.
 * TODO: more sophisticated check for bond structure (e.g. connected subgraph, ...)
 * \return true - bonds present, false - no bonds
 */
bool molecule::hasBondStructure() const
{
  for(molecule::const_iterator AtomRunner = begin(); AtomRunner != end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    if (!ListOfBonds.empty())
      return true;
  }
  return false;
}

/** Prints a list of all bonds to \a *out.
 */
void molecule::OutputBondsList() const
{
  DoLog(1) && (Log() << Verbose(1) << endl << "From contents of bond chain list:");
  for(molecule::const_iterator AtomRunner = molecule::begin(); AtomRunner != molecule::end(); ++AtomRunner) {
    const BondList& ListOfBonds = (*AtomRunner)->getListOfBonds();
    for(BondList::const_iterator BondRunner = ListOfBonds.begin();
        BondRunner != ListOfBonds.end();
        ++BondRunner)
      if ((*BondRunner)->leftatom == *AtomRunner) {
        DoLog(0) && (Log() << Verbose(0) << *(*BondRunner) << "\t" << endl);
      }
  }
  DoLog(0) && (Log() << Verbose(0) << endl);
}
;


/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::Nr and then the Atom::Nr of all bond partners per line.
 * \param &filename name of file
 * \param path path to file, defaults to empty
 * \return true - file written successfully, false - writing failed
 */
bool molecule::StoreAdjacencyToFile(std::string filename, std::string path)
{
  ofstream AdjacencyFile;
  string line;
  bool status = true;

  if (path != "")
    line = path + "/" + filename;
  else
    line = filename;
  AdjacencyFile.open(line.c_str(), ios::out);
  DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
  if (AdjacencyFile.good()) {
    AdjacencyFile << "m\tn" << endl;
    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputAdjacency),&AdjacencyFile));
    AdjacencyFile.close();
    DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
  } else {
    DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
    status = false;
  }

  return status;
}
;

/** Storing the bond structure of a molecule to file.
 * Simply stores Atom::Nr and then the Atom::Nr of all bond partners, one per line.
 * \param &filename name of file
 * \param path path to file, defaults to empty
 * \return true - file written successfully, false - writing failed
 */
bool molecule::StoreBondsToFile(std::string filename, std::string path)
{
  ofstream BondFile;
  string line;
  bool status = true;

  if (path != "")
    line = path + "/" + filename;
  else
    line = filename;
  BondFile.open(line.c_str(), ios::out);
  DoLog(1) && (Log() << Verbose(1) << "Saving adjacency list ... " << endl);
  if (BondFile.good()) {
    BondFile << "m\tn" << endl;
    for_each(atoms.begin(),atoms.end(),bind2nd(mem_fun(&atom::OutputBonds),&BondFile));
    BondFile.close();
    DoLog(1) && (Log() << Verbose(1) << "\t... done." << endl);
  } else {
    DoLog(1) && (Log() << Verbose(1) << "\t... failed to open file " << line << "." << endl);
    status = false;
  }

  return status;
}
;

bool CheckAdjacencyFileAgainstMolecule_Init(std::string &path, ifstream &File, int *&CurrentBonds)
{
  string filename;
  filename = path + ADJACENCYFILE;
  File.open(filename.c_str(), ios::out);
  DoLog(1) && (Log() << Verbose(1) << "Looking at bond structure stored in adjacency file and comparing to present one ... " << endl);
  if (File.fail())
    return false;

  // allocate storage structure
  CurrentBonds = new int[MAXBONDS]; // contains parsed bonds of current atom
  for(int i=0;i<MAXBONDS;i++)
    CurrentBonds[i] = 0;
  return true;
}
;

void CheckAdjacencyFileAgainstMolecule_Finalize(ifstream &File, int *&CurrentBonds)
{
  File.close();
  File.clear();
  delete[](CurrentBonds);
}
;

void CheckAdjacencyFileAgainstMolecule_CompareBonds(bool &status, int &NonMatchNumber, atom *&Walker, size_t &CurrentBondsOfAtom, int AtomNr, int *&CurrentBonds, atom **ListOfAtoms)
{
  size_t j = 0;
  int id = -1;

  //Log() << Verbose(2) << "Walker is " << *Walker << ", bond partners: ";
  const BondList& ListOfBonds = Walker->getListOfBonds();
  if (CurrentBondsOfAtom == ListOfBonds.size()) {
    for (BondList::const_iterator Runner = ListOfBonds.begin();
        Runner != ListOfBonds.end();
        ++Runner) {
      id = (*Runner)->GetOtherAtom(Walker)->getNr();
      j = 0;
      for (; (j < CurrentBondsOfAtom) && (CurrentBonds[j++] != id);)
        ; // check against all parsed bonds
      if (CurrentBonds[j - 1] != id) { // no match ? Then mark in ListOfAtoms
        ListOfAtoms[AtomNr] = NULL;
        NonMatchNumber++;
        status = false;
        DoeLog(2) && (eLog() << Verbose(2) << id << " can not be found in list." << endl);
      } else {
        //Log() << Verbose(0) << "[" << id << "]\t";
      }
    }
    //Log() << Verbose(0) << endl;
  } else {
    DoLog(0) && (Log() << Verbose(0) << "Number of bonds for Atom " << *Walker << " does not match, parsed " << CurrentBondsOfAtom << " against " << ListOfBonds.size() << "." << endl);
    status = false;
  }
}
;

/** Checks contents of adjacency file against bond structure in structure molecule.
 * \param *path path to file
 * \param **ListOfAtoms allocated (molecule::AtomCount) and filled lookup table for ids (Atom::Nr) to *Atom
 * \return true - structure is equal, false - not equivalence
 */
bool molecule::CheckAdjacencyFileAgainstMolecule(std::string &path, atom **ListOfAtoms)
{
  ifstream File;
  bool status = true;
  atom *Walker = NULL;
  int *CurrentBonds = NULL;
  int NonMatchNumber = 0; // will number of atoms with differing bond structure
  size_t CurrentBondsOfAtom = -1;
  const int AtomCount = getAtomCount();

  if (!CheckAdjacencyFileAgainstMolecule_Init(path, File, CurrentBonds)) {
    DoLog(1) && (Log() << Verbose(1) << "Adjacency file not found." << endl);
    return true;
  }

  char buffer[MAXSTRINGSIZE];
  int tmp;
  // Parse the file line by line and count the bonds
  while (!File.eof()) {
    File.getline(buffer, MAXSTRINGSIZE);
    stringstream line;
    line.str(buffer);
    int AtomNr = -1;
    line >> AtomNr;
    CurrentBondsOfAtom = -1; // we count one too far due to line end
    // parse into structure
    if ((AtomNr >= 0) && (AtomNr < AtomCount)) {
      Walker = ListOfAtoms[AtomNr];
      while (line >> ws >> tmp) {
        std::cout << "Recognized bond partner " << tmp << std::endl;
        CurrentBonds[++CurrentBondsOfAtom] = tmp;
        ASSERT(CurrentBondsOfAtom < MAXBONDS,
            "molecule::CheckAdjacencyFileAgainstMolecule() - encountered more bonds than allowed: "
            +toString(CurrentBondsOfAtom)+" >= "+toString(MAXBONDS)+"!");
      }
      // compare against present bonds
      CheckAdjacencyFileAgainstMolecule_CompareBonds(status, NonMatchNumber, Walker, CurrentBondsOfAtom, AtomNr, CurrentBonds, ListOfAtoms);
    } else {
      if (AtomNr != -1)
        DoeLog(2) && (eLog() << Verbose(2) << AtomNr << " is not valid in the range of ids [" << 0 << "," << AtomCount << ")." << endl);
    }
  }
  CheckAdjacencyFileAgainstMolecule_Finalize(File, CurrentBonds);

  if (status) { // if equal we parse the KeySetFile
    DoLog(1) && (Log() << Verbose(1) << "done: Equal." << endl);
  } else
    DoLog(1) && (Log() << Verbose(1) << "done: Not equal by " << NonMatchNumber << " atoms." << endl);
  return status;
}
;

/** Adds a bond as a copy to a given one
 * \param *left leftatom of new bond
 * \param *right rightatom of new bond
 * \param *CopyBond rest of fields in bond are copied from this
 * \return pointer to new bond
 */
bond * molecule::CopyBond(atom *left, atom *right, bond *CopyBond)
{
  bond *Binder = AddBond(left, right, CopyBond->BondDegree);
  Binder->Cyclic = CopyBond->Cyclic;
  Binder->Type = CopyBond->Type;
  return Binder;
}
;

void BuildInducedSubgraph_Init(atom **&ParentList, int AtomCount)
{
  // reset parent list
  ParentList = new atom*[AtomCount];
  for (int i=0;i<AtomCount;i++)
    ParentList[i] = NULL;
  DoLog(3) && (Log() << Verbose(3) << "Resetting ParentList." << endl);
}
;

void BuildInducedSubgraph_FillParentList(const molecule *mol, const molecule *Father, atom **&ParentList)
{
  // fill parent list with sons
  DoLog(3) && (Log() << Verbose(3) << "Filling Parent List." << endl);
  for (molecule::const_iterator iter = mol->begin(); iter != mol->end(); ++iter) {
    ParentList[(*iter)->father->getNr()] = (*iter);
    // Outputting List for debugging
    DoLog(4) && (Log() << Verbose(4) << "Son[" << (*iter)->father->getNr() << "] of " << (*iter)->father << " is " << ParentList[(*iter)->father->getNr()] << "." << endl);
  }
};

void BuildInducedSubgraph_Finalize(atom **&ParentList)
{
  delete[](ParentList);
}
;

bool BuildInducedSubgraph_CreateBondsFromParent(molecule *mol, const molecule *Father, atom **&ParentList)
{
  bool status = true;
  atom *OtherAtom = NULL;
  // check each entry of parent list and if ok (one-to-and-onto matching) create bonds
  DoLog(3) && (Log() << Verbose(3) << "Creating bonds." << endl);
  for (molecule::const_iterator iter = Father->begin(); iter != Father->end(); ++iter) {
    if (ParentList[(*iter)->getNr()] != NULL) {
      if (ParentList[(*iter)->getNr()]->father != (*iter)) {
        status = false;
      } else {
        const BondList& ListOfBonds = (*iter)->getListOfBonds();
        for (BondList::const_iterator Runner = ListOfBonds.begin();
            Runner != ListOfBonds.end();
            ++Runner) {
          OtherAtom = (*Runner)->GetOtherAtom((*iter));
          if (ParentList[OtherAtom->getNr()] != NULL) { // if otheratom is also a father of an atom on this molecule, create the bond
            DoLog(4) && (Log() << Verbose(4) << "Endpoints of Bond " << (*Runner) << " are both present: " << ParentList[(*iter)->getNr()]->getName() << " and " << ParentList[OtherAtom->getNr()]->getName() << "." << endl);
            mol->AddBond(ParentList[(*iter)->getNr()], ParentList[OtherAtom->getNr()], (*Runner)->BondDegree);
          }
        }
      }
    }
  }
  return status;
}
;

/** Adds bond structure to this molecule from \a Father molecule.
 * This basically causes this molecule to become an induced subgraph of the \a Father, i.e. for every bond in Father
 * with end points present in this molecule, bond is created in this molecule.
 * Special care was taken to ensure that this is of complexity O(N), where N is the \a Father's molecule::AtomCount.
 * \param *Father father molecule
 * \return true - is induced subgraph, false - there are atoms with fathers not in \a Father
 * \todo not checked, not fully working probably
 */
bool molecule::BuildInducedSubgraph(const molecule *Father){
  bool status = true;
  atom **ParentList = NULL;
  DoLog(2) && (Log() << Verbose(2) << "Begin of BuildInducedSubgraph." << endl);
  BuildInducedSubgraph_Init(ParentList, Father->getAtomCount());
  BuildInducedSubgraph_FillParentList(this, Father, ParentList);
  status = BuildInducedSubgraph_CreateBondsFromParent(this, Father, ParentList);
  BuildInducedSubgraph_Finalize(ParentList);
  DoLog(2) && (Log() << Verbose(2) << "End of BuildInducedSubgraph." << endl);
  return status;
}
;

/** Fills a lookup list of father's Atom::nr -> atom for each subgraph.
 * \param **&ListOfLocalAtoms Lookup table for each subgraph and index of each atom in global molecule, may be NULL on start, then it is filled
 * \param GlobalAtomCount number of atoms in the complete molecule
 * \return true - success, false - failure (ListOfLocalAtoms != NULL)
 */
bool molecule::FillListOfLocalAtoms(atom **&ListOfLocalAtoms, const int GlobalAtomCount)
{
  bool status = true;

  if (ListOfLocalAtoms == NULL) { // allocate and fill list of this fragment/subgraph
    status = status && CreateFatherLookupTable(ListOfLocalAtoms, GlobalAtomCount);
  } else
    return false;

  return status;
}