timvdm

#include <iostream>
#include <vector>
#include <algorithm>
#include <cassert>

#include <Eigen/Core>

#include <openbabel/mol.h>
#include <openbabel/obconversion.h>
#include <openbabel/graphsym.h>

#include <graph.hh>

#include "permutation.h"
using namespace OpenBabel;

////////////////////////////////////////////////////////////////////////////////
//
//
//  O P E N B A B E L   S P E C I F I C   C O D E
//
//
////////////////////////////////////////////////////////////////////////////////

OpenBabel::OBMol* ReadMolecule(const std::string &filename)
{
  OpenBabel::OBConversion conv;
  OpenBabel::OBFormat *format = conv.FormatFromExt(filename.c_str());
  if (!format || !conv.SetInFormat(format)) {
    std::cout << "ERROR: Could not detect file format for: " << filename << std::endl;
    return 0;
  }

  std::ifstream ifs;
  ifs.open(filename.c_str());
  if (!ifs) {
    std::cout << "ERROR: Could not open file: " << filename << std::endl;
    return 0;
  }

  OpenBabel::OBMol *obmol = new OpenBabel::OBMol;
  if (!conv.Read(obmol, &ifs)) {
    std::cout << "ERROR: Could not read molecule from file: " << filename << std::endl;
    delete obmol;
    ifs.close();
    return 0;
  }

  ifs.close();
  return obmol;
}

////////////////////////////////////////////////////////////////////////////////
//
//
//  U N I T   T E S T S
//
//
////////////////////////////////////////////////////////////////////////////////



void testPermutation()
{
  std::vector<unsigned int> e;
  e.push_back(1);
  e.push_back(2);
  e.push_back(3);
  e.push_back(4);
  Permutation p1(e), p2(e);

  assert(p1 == p2);
}

void testPermutationGroup()
{
  Permutation p;
  p.map.push_back(1);
  p.map.push_back(2);
  p.map.push_back(3);
  p.map.push_back(4);

  PermutationGroup group;
  assert( !group.contains(p) );

  group.add(p);
  assert( group.contains(p) );
}

void unitTests()
{
  testPermutation();
  testPermutationGroup();

}




void callback(void *param, unsigned int n, const unsigned int *aut)
{
  assert(param);
  /*
  fprintf((FILE*)param, "Generator: ");
  bliss::print_permutation((FILE*)param, n, aut, 0);
  fprintf((FILE*)param, "\n");
  */

  std::vector<unsigned int> elements;
  cout << "Generator: ";
  for (unsigned int i = 0; i < n; ++i) {
    cout << aut[i] + 1 << " ";
    elements.push_back(aut[i]+1);
  }
  cout << endl;

  PermutationGroup *generators = static_cast<PermutationGroup*>(param);
  generators->add(Permutation(elements));
}

void addInverses(PermutationGroup &G)
{
  unsigned int size = G.permutations.size();
  for (unsigned int i = 0; i < size; ++i) {
    Permutation inv_p(G.permutations.at(i).matrix().transpose());
    if (!G.contains(inv_p))
      G.add(inv_p);
  }
}

void addProducts(PermutationGroup &G)
{
  for (unsigned int i = 0; i < G.permutations.size(); ++i) {
    for (unsigned int j = 0; j < G.permutations.size(); ++j) {
      if (i >= j)
        continue;

      Permutation p = G.permutations.at(i) * G.permutations.at(j);
      if (!G.contains(p))
        G.add(p);
    }
  }
}

PermutationGroup findAutomorphisms(OpenBabel::OBMol *obmol, const std::vector<unsigned int> &symClasses)
{
  using OpenBabel::OBMolAtomIter;
  using OpenBabel::OBMolBondIter;

  // construct the bliss graph
  bliss::Graph g;
  std::map<OpenBabel::OBAtom*, unsigned int> atom2index;
  FOR_ATOMS_OF_MOL (atom, obmol) {
    atom2index[&*atom] = g.add_vertex(symClasses.at(atom->GetIndex()));
  }
  FOR_BONDS_OF_MOL (bond, obmol) {
    g.add_edge(atom2index[bond->GetBeginAtom()], atom2index[bond->GetEndAtom()]);
  }

  // use bliss to get the automorphism group generators
  PermutationGroup generators;
  bliss::Stats stats;
  g.find_automorphisms(stats, &callback, &generators);

  cout << "# Automorphisms:" << stats.group_size_approx << endl;
  unsigned long nAut = stats.group_size_approx;

  // construct the automorphism group
  PermutationGroup G;

  // add identity permutation
  std::vector<unsigned int> eElements;
  for (unsigned int i = 0; i < obmol->NumAtoms(); ++i)
    eElements.push_back(i+1);
  Permutation e(eElements);
  if (!G.contains(e))
    G.add(e);

  // add the generators
  for (unsigned int i = 0; i < generators.size(); ++i) {
    if (!G.contains(generators.at(i)))
      G.add(generators.at(i));
  }


  unsigned int counter = 0;
  unsigned int lastSize = 0;
  while (G.permutations.size() != lastSize) {
    lastSize = G.permutations.size();

    addInverses(G);
    addProducts(G);

    counter++;
    if (counter > 100)
      break;
  }

  return G;
}






int main(int argc, char **argv)
{
  unitTests();

  if (argc < 2) {
    std::cout << "Usage: " << argv[0] << " <filename>" << std::endl;
    return -1;
  }

  OpenBabel::OBMol *obmol = ReadMolecule(argv[1]);
  if (!obmol)
    return -1;

  std::vector<unsigned int> elements;
  for (int i = 0; i < obmol->NumAtoms(); ++i) {
    elements.push_back(i+1);
  }

  OpenBabel::OBGraphSym graphSym(obmol);
  std::vector<unsigned int> symClasses;
  graphSym.GetSymmetry(symClasses);


  PermutationGroup G = findAutomorphisms(obmol, symClasses);

  std::cout << "Automorphisms:" << std::endl;
  for (unsigned int i = 0; i < G.size(); ++i)
    G.at(i).print();


  return 0;
}