timvdm

/***************************************************************************
 *   Copyright (C) 2009 by Tim Vandermeersch                               *
 *                                                                         *
 *   This program 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.                                   *
 *                                                                         *
 *   This program 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 this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include <iostream>
#include <vector>
#include <algorithm>
#include <cassert>

#include <Eigen/Core>

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

/**
 * This class represents a shortened mapping for permutations.
 *
 * There are three ways to symbolically represent permutations:
 @verbatim
 cycle notation:  (1 2 3)    (rotate to the left: 1 2 3 -> 2 3 1)

 standard mapping:  1 2 3    (first line are the original numbers in ascending order)
                    2 3 1    (second line contains the element to which it is permuted)

 shortened mapping:  2 3 1   (avoid first line in standard mapping: efficient storage)
 @endverbatim
 *
 * For efficiency reasons, the shortened mapping notation is used to internally store
 * the permutations. The elements are always in the range \f$0-n\f$ where \f$n\f$ is the
 * number of elements (i.e. atoms).
 *
 * Multiple Permutations can be combined in a PermutationGroup. The permutation
 * group \f$S_n\f$, is the group containing all \f$n!\f$ permutations.
 *
 * The automorphism group is a subgroup of \f$S_n\f$. It consists of all vertex (i.e. atom)
 * permutations which conserve the edge connections (i.e. bonds). There is always
 * at least one permutation in the automorphism group, namely the identity permutation.
 * Since it doesn't change any connections, it fulfills the previously given definition.
 * To check if other permuations belong to the automorphism group, the following matrix
 * multiplication equation can be used
 \f[
   P^T A P = A
 \f]
 * For small molecules (less than 10 atoms) it is possible to generate all the
 * permutations in the \f$S_n\f$ group and apply the formula. However, to have reasonable
 * performance for larger molecules, it is needed to limit the number of permutations
 * that need to be checked. In practise, a Graph Invariant Atom Partitioning (GIAP)
 * algorithm is used for this (e.g. Morgan's Extended Connectivity (EC)). See the
 * OBGraphSym class documentation for more details. [note: an orbit is the group of atoms
 * in a molecule that belong to the same partitioned group -- or the atoms with the same
 * symmetry class]
 *
 * Once the initial atom partitioning in orbits is obtained, this can be used to generate
 * the subgroup of \f$S_n\f$ knowing that permutations between orbits are not allowed
 * (i.e. these permutations will never satisfy the matrix multiplication equation).
 * Orbits with only one atom will never be changed in any permutation belonging to the
 * automorphisms group and will therefor be ignored. The number of permutations to be
 * tested depends on the number of remaining orbits and the number of particles contained
 * in them. For example a molecule with 12 atoms, 2 orbits, each containing 3 atoms will require
 * \f$3! \times 3! = 36\f$ tests. [note: \f$12! = 4790 001 600\f$]
 *
 * @image html 12a2o3a.png
 *
 */
struct Permutation
{
  /**
   * Default constructor.
   */
  Permutation()
  {}

  /**
   * Constructor taking a @p map as argument.
   */
  Permutation(const std::vector<unsigned int> &_map) : map(_map)
  {}

  /**
   * Constructor taking a @p matrix as argument.
   */
  Permutation(const Eigen::MatrixXi &matrix)
  {
    setMatrix(matrix);
  }

  /**
   * Copy constructor.
   */
  Permutation(const Permutation &other)
  {
    this->map = other.map;
  }

  /**
   * Print the permutation to std::cout in shortened notation.
   */
  void print() const
  {
    std::vector<unsigned int>::const_iterator i;
    for (i = map.begin(); i != map.end(); ++i)
      std::cout << *i << " ";
    std::cout << std::endl;
  }

  /**
   * Compute the permutation matrix for this Permutation.
   */
  Eigen::MatrixXi matrix() const
  {
    int n = map.size();
    Eigen::MatrixXi P = Eigen::MatrixXi::Zero(n, n);

    for (int i = 0; i < n; ++i) {
      P(map.at(i)-1, i) = 1;
      P(i, map.at(i)-1) = 1;
    }

    return P;
  }

  /**
   * Set the permutation matrix for this permutation. The matrix is
   * converted to a shortened notation permutation.
   */
  void setMatrix(const Eigen::MatrixXi &m)
  {
    if (m.rows() != m.cols())
      return;

    int size = m.rows();
    map.resize(size);

    for (int i = 0; i < size; ++i)
      for (int j = 0; j < size; ++j)
        if (m(i,j))
          map[i] = j + 1;
  }

  /**
   * Multiply two permutation matrices and return the resulting permutation.
   * This is the same as applying the two permutations consecutively.
   */
  Permutation operator*(const Permutation &rhs) const
  {
    Eigen::MatrixXi res = rhs.matrix() * matrix();
    return Permutation(res);
  }

  /**
   * Equality operator. Two permutations are equal if their shortened notations
   * are the same (compared element-by-element).
   */
  bool operator==(const Permutation &rhs) const
  {
    if (map.size() != rhs.map.size())
      return false;

    std::vector<unsigned int>::const_iterator i1 = map.begin();
    std::vector<unsigned int>::const_iterator i2 = rhs.map.begin();
    for (; i1 != map.end(); ++i1, ++i2)
      if (*i1 != *i2)
        return false;

    return true;
  }

  std::vector<unsigned int> map; //!< The actual mapping in shortened notation
};

struct PermutationGroup
{
  /**
   * Default constructor.
   */
  PermutationGroup()
  {}

  /**
   * Constructor taking vector of permutations as argument.
   */
  PermutationGroup(const std::vector<Permutation> &_permutations) : permutations(_permutations)
  {}

  /**
   * @return The number of permutations in this group.
   */
  unsigned int size() const
  {
    return permutations.size();
  }

  /**
   * Add permutation @p p to this permutation group.
   */
  void add(const Permutation &p)
  {
    permutations.push_back(p);
  }

  /**
   * @return A constant reference to the @p index-th permutation.
   */
  const Permutation& at(unsigned int index)
  {
    return permutations.at(index);
  }

  /**
   * @return True if this permutation group contains permutation @p p.
   */
  bool contains(const Permutation &p) const
  {
    std::vector<Permutation>::const_iterator i;
    for (i = permutations.begin(); i != permutations.end(); ++i)
      if (*i == p)
        return true;
    return false;
  }

  std::vector<Permutation> permutations; //!< The actual permutations in the group
};


struct Orbit
{
  Orbit(const std::vector<unsigned int> &_ordered) : ordered(_ordered), current(ordered)
  {
  }

  Orbit(const Orbit &other) : ordered(other.ordered), current(other.ordered)
  {
  }

  std::vector<unsigned int> ordered;
  std::vector<unsigned int> current;
};

typedef std::vector<Orbit> Orbits;

/**
 * Compute the \f$S_n\f$ group consisting of all \f$n!\f$ permutations.
 */
PermutationGroup computeAllPermutations(const Permutation &input)
{
  PermutationGroup Sn;
  Permutation p = input;
  Sn.add(p);

  while (std::next_permutation(p.map.begin(), p.map.end())) {
    Sn.add(p);
  }

  return Sn;
}

Orbits computeOrbits(const std::vector<unsigned int> &symClasses)
{
  Orbits orbits;
  std::vector<bool> visit(symClasses.size(), false);

  std::vector<unsigned int>::const_iterator symClass;
  for (symClass = symClasses.begin(); symClass != symClasses.end(); ++symClass) {
    if (visit.at(*symClass))
      continue;
    visit[*symClass] = true;

    unsigned int n = std::count(symClasses.begin(), symClasses.end(), *symClass);
    if (n >= 2) {
      std::vector<unsigned int> orbit;

      for (unsigned int i = 0; i < symClasses.size(); ++i)
        if (symClasses.at(i) == *symClass)
          orbit.push_back(i);

      orbits.push_back(Orbit(orbit));
    }
  }

  return orbits;
}

Permutation processOrbitPermutations(const Orbit &orbit, const Permutation &E)
{
  Permutation p = E;
  for (unsigned int i = 0; i < orbit.ordered.size(); ++i) {
    std::cout << "p.map.at(" << orbit.ordered.at(i) << ") = E.map.at(" << orbit.current.at(i) << ")" << std::endl;
    p.map.at(orbit.ordered.at(i)) = E.map.at(orbit.current.at(i));
  }
  return p;
}

bool sameMatrix(const Eigen::MatrixXi &m1, const Eigen::MatrixXi &m2)
{
  if (m1.cols() != m2.cols())
    return false;
  if (m1.rows() != m2.rows())
    return false;

  for (int i = 0; i < m1.rows(); ++i)
    for (int j = 0; j < m1.cols(); ++j)
      if (m1(j,i) != m2(j,i))
        return false;

  return true;
}

PermutationGroup computeAutomorphisms(const PermutationGroup &g, const Eigen::MatrixXi &A)
{
  PermutationGroup G;

  std::vector<Permutation>::const_iterator p;
  for (p = g.permutations.begin(); p != g.permutations.end(); ++p) {
    Eigen::MatrixXi P = (*p).matrix();
    Eigen::MatrixXi A2 = P.transpose() * A *  P;

    if (sameMatrix(A, A2)) {
      G.add(*p);
      (*p).print();
    }
  }

  return G;
}


/**
 * Compute the adjacency matrix for the specified molecule. This matrix
 * represents the connectivity in the molecule. It can be used to determine
 * if a Permutation is part of the automorphism group.
 */
Eigen::MatrixXi adjacencyMatrix(OpenBabel::OBMol *mol)
{
  int n = mol->NumAtoms();
  Eigen::MatrixXi A = Eigen::MatrixXi::Zero(n, n);

  using OpenBabel::OBMolBondIter;
  std::cout << "bonds:" << std::endl;
  FOR_BONDS_OF_MOL (bond, mol)
    std::cout << "  " << bond->GetBeginAtomIdx() << "-" << bond->GetEndAtomIdx() << std::endl;

  for (int i = 0; i < n; ++i) {
    OpenBabel::OBAtom *ai = mol->GetAtom(i+1);
    for (int j = 0; j < n; ++j) {
      if (i == j)
        continue;
      if (ai->IsConnected(mol->GetAtom(j+1)))
        A(i, j) = 1;
    }
  }

  return A;
}

////////////////////////////////////////////////////////////////////////////////
//
//
//  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();

}

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);
  }

  Permutation E(elements); // identity

  OpenBabel::OBGraphSym graphSym(obmol);
  std::vector<unsigned int> symClasses;
  graphSym.GetSymmetry(symClasses);
  const Orbits sortedOrbits = computeOrbits(symClasses);

  PermutationGroup G0, G1;

  Orbits orbits = sortedOrbits;;
  for (unsigned int i = 0; i < orbits.size(); ++i) {
    // print the orbit
    std::cout << "orbit: ";
    for (unsigned int j = 0; j < orbits.at(i).ordered.size(); ++j) {
      std::cout << orbits.at(i).ordered.at(j) << " ";
    }
    std::cout << std::endl;

    // sorted orbit
    Orbit e = orbits.at(i);

    Permutation p = processOrbitPermutations(orbits.at(i), E);
    if (!G0.contains(p))
      G0.add(p);
    while (std::next_permutation(orbits.at(i).current.begin(), orbits.at(i).current.end())) {
      p = processOrbitPermutations(orbits.at(i), E);

      if (!G0.contains(p))
        G0.add(p);
    }
    std::cout << std::endl;
  }

  for (unsigned int i = 0; i < G0.size(); ++i) {
    for (unsigned int j = 0; j < G0.size(); ++j) {
      if (i == j)
        continue;

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

  for (unsigned int i = 0; i < G1.size(); ++i)
    if (!G0.contains(G1.at(i)))
      G0.add(G1.at(i));

  G1.permutations.clear();

  // print automorphisms
  std::cout << "G0:" << std::endl;
  for (unsigned int i = 0; i < G0.size(); ++i)
    G0.at(i).print();
  // print automorphisms
  std::cout << "G1:" << std::endl;
  for (unsigned int i = 0; i < G1.size(); ++i)
    G1.at(i).print();


  // compute adjacency matrix
  Eigen::MatrixXi A = adjacencyMatrix(obmol);
  // compute the permutation group of all n! permutations
  PermutationGroup Sn = computeAllPermutations(E);
  // compute automorphisms starting with the limited group
  PermutationGroup G3 = computeAutomorphisms(G0, A);
  // compute automorphisms starting with entire permutation group
  PermutationGroup G4 = computeAutomorphisms(Sn, A);


  if (G3.size() != G4.size()) {
    std::cout << "NOT ALL AUTOMORPHISMS ARE FOUND!!" << std::endl;
    return 0;
  }
  bool success = true;
  for (unsigned int i = 0; i < G4.size(); ++i)
    if (!G3.contains(G4.at(i))) {
      std::cout << "COULD NOT FIND AUTOMORPHISM: ";
      G3.at(i).print();
      std::cout << std::endl;
      success = false;
    }

  if (success) {
    std::cout << "ALL AUTOMORPHISMS FOUND!!" << std::endl;
  }

  return 0;
}