// $Id$
//
// Copyright (C) 2013-2014 Paolo Tosco
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#include "O3AAlignMolecules.h"
#include <RDGeneral/utils.h>
#include <Numerics/Vector.h>
#include <GraphMol/Substruct/SubstructMatch.h>
#include <GraphMol/Conformer.h>
#include <GraphMol/ROMol.h>
#include <GraphMol/MolOps.h>
#include <GraphMol/AtomIterators.h>
#include <Numerics/Alignment/AlignPoints.h>
#include <GraphMol/MolTransforms/MolTransforms.h>
#include <boost/dynamic_bitset.hpp>
#include <boost/algorithm/string.hpp>
#include <boost/math/special_functions/round.hpp>

#define square(x) ((x) * (x))


namespace RDKit {
  namespace MolAlign {
    static boost::uint8_t mmffSimMatrix[99][99] = {
      {  1,  3,  4,  3,  0,  4,  6,  7,  5,  6,  7,  4,  5,  7,  3,  5,  6,  7,  3,  2,  0,  2,  0,  0,  6,  7,  0,  0,  0,  3,  0, 10,  0, 10, 10,  0,  3,  6,  7,  6,  8,  7,  8,  4,  7,  6,  8,  7, 10,  0, 10,  0,  6, 10, 12, 14, 12, 10,  4,  4,  6, 10,  3,  3,  6,  6,  8,  8,  8,  8,  0, 10, 10,  4,  5, 10, 10,  3,  6, 10, 12,  8,  0,  0,  0,  0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  3,  1,  2,  2,  0,  5,  6,  7,  3,  5,  7,  4,  5,  7,  3,  4,  5,  6,  3,  2,  0,  2,  0,  0,  5,  8,  0,  0,  0,  2,  0, 10,  0,  8, 10,  0,  2,  4,  5,  5,  6,  4,  7,  4,  6,  5,  7,  6, 10,  0,  8,  0,  4, 10, 12, 14, 12, 10,  4,  3,  5, 10,  2,  2,  4,  4,  7,  7,  7,  8,  0, 10,  8,  3,  4, 10, 10,  2,  4, 10, 12,  7,  0,  0,  0,  0,  8,  8, 10, 10, 10,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  4,  2,  1,  2,  0,  7,  7,  8,  2,  5,  7,  4,  5,  7,  5,  5,  3,  4,  4,  4,  0,  4,  0,  0,  4,  8,  0,  0,  0,  2,  0, 10,  0,  8, 12,  0,  2,  4,  6,  6,  4,  8,  8,  6,  6,  3,  8,  7,  8,  0,  6,  0,  6,  8, 10, 12, 10,  8,  6,  6,  3, 12,  2,  2,  6,  6,  4,  4,  4, 10,  0, 12,  5,  3,  6, 12, 12,  2,  4,  7, 10,  5,  0,  0,  0,  0,  6,  6, 12, 12, 12,  6,  6,  6,  6,  6,  6,  6,  6 },
      {  3,  2,  2,  1,  0,  5,  6,  7,  4,  5,  7,  4,  5,  7,  3,  4,  5,  6,  3,  2,  0,  2,  0,  0,  5,  8,  0,  0,  0,  2,  0, 10,  0,  8, 10,  0,  2,  4,  5,  5,  6,  4,  7,  4,  6,  5,  7,  6, 10,  0,  8,  0,  4, 10, 12, 14, 12, 10,  4,  3,  5, 10,  2,  2,  4,  4,  7,  7,  7,  8,  0, 10,  8,  3,  4, 10, 10,  2,  4, 10, 12,  7,  0,  0,  0,  0,  8,  8, 10, 10, 10,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  4,  5,  7,  5,  0,  1,  3,  5,  4,  5,  6,  5,  5,  6,  2,  4,  7,  7,  5,  4,  0,  4,  0,  0,  7,  6,  0,  0,  0,  4,  0,  5,  0, 12,  5,  0,  6,  6,  4,  4,  6,  4,  4,  4,  7,  6,  6,  6,  4,  0,  6,  0,  5, 12, 14, 16, 14, 12,  2,  4,  8,  6,  6,  6,  4,  4,  8,  8,  8,  3,  0,  6,  6,  7,  6,  6,  8,  5,  4, 12, 10,  8,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  6,  7,  6,  0,  3,  1,  6,  3,  5,  6,  5,  5,  6,  4,  2,  7,  7,  7,  6,  0,  6,  0,  0,  7,  6,  0,  0,  0,  4,  0,  4,  0, 12,  4,  0,  6,  6,  4,  4,  7,  4,  4,  5, 10,  8,  6,  5, 10,  0, 10,  0,  5, 12, 14, 16, 14, 12,  3,  4,  8,  6,  6,  6,  4,  4,  8,  8,  8,  4,  0,  6,  6,  8,  5,  6,  8,  6,  6, 12, 10,  8,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  7,  7,  8,  7,  0,  5,  6,  1,  2,  6,  8,  7,  7,  8,  6,  7,  7,  7,  7,  7,  0,  7,  0,  0,  7,  5,  0,  0,  0,  7,  0, 12,  0,  4, 12,  0,  7,  2,  4,  2, 12,  3,  6,  7,  8,  8,  7,  7, 10,  0, 10,  0,  5,  4,  5,  5,  8,  4,  6,  6,  6,  5,  7,  7,  4,  4,  7,  7,  7,  6,  0, 10, 10, 10,  6,  6, 10,  7,  4,  7,  5,  7,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  5,  3,  2,  4,  0,  4,  3,  2,  1,  5,  8,  7,  7,  8,  6,  6,  7,  7,  7,  5,  0,  5,  0,  0,  7,  6,  0,  0,  0,  5,  0, 10,  0,  5, 10,  0,  5,  2,  4,  2, 10,  2,  5,  6,  8,  7,  5,  4, 10,  0,  8,  0,  4,  3,  3,  3,  6,  3,  5,  5,  5,  4,  5,  5,  4,  4,  7,  7,  7,  7,  0, 10, 10, 10,  3,  7, 10,  5,  4,  6,  4,  7,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  6,  5,  5,  5,  0,  5,  5,  6,  5,  1,  8,  7,  7,  8,  5,  7,  8,  8,  8,  6,  0,  6,  0,  0,  8,  6,  0,  0,  0,  6,  0,  6,  0,  8,  6,  0,  6,  3,  3,  5,  8,  3,  2,  6,  8,  7,  4,  5, 12,  0, 10,  0,  5,  6,  6,  8, 10,  5,  6,  6,  6,  3,  6,  6,  4,  4,  7,  7,  7,  8,  0,  8,  8, 10,  6,  4, 10,  6,  4,  8,  6,  7,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  7,  7,  7,  7,  0,  6,  6,  8,  8,  8,  1,  3,  4,  6,  6,  6, 10, 10, 10,  8,  0,  8,  0,  0,  8, 10,  0,  0,  0,  8,  0,  3,  0, 12,  3,  0,  7,  7,  7,  8,  4,  5,  6,  6, 10,  8,  8,  6, 12,  0, 12,  0, 10, 12, 12, 14, 12, 12,  5,  6,  8,  5,  7,  7,  6,  6, 10, 10, 10,  5,  0,  4,  6, 10,  6,  3,  5,  7,  6, 12, 14, 10,  0,  0,  0,  0, 14, 14,  3,  4,  5, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  4,  4,  4,  4,  0,  5,  5,  7,  7,  7,  3,  1,  2,  3,  4,  4,  8,  8,  6,  6,  0,  6,  0,  0,  7,  8,  0,  0,  0,  6,  0,  4,  0, 10,  4,  0,  5,  5,  5,  6,  5,  4,  5,  4,  8,  7,  7,  6, 10,  0, 10,  0, 10, 12, 12, 14, 12, 12,  4,  5,  8,  6,  5,  5,  5,  5, 10, 10, 10,  6,  0,  4,  6, 10,  6,  4,  4,  5,  5, 12, 14, 10,  0,  0,  0,  0, 12, 12,  4,  3,  3, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  5,  5,  5,  5,  0,  5,  5,  7,  7,  7,  4,  2,  1,  2,  4,  4,  8,  8,  6,  6,  0,  6,  0,  0,  7,  8,  0,  0,  0,  6,  0,  4,  0, 10,  4,  0,  5,  5,  5,  6,  5,  4,  5,  4,  8,  7,  7,  6, 10,  0, 10,  0, 10, 12, 12, 14, 12, 12,  4,  5,  8,  6,  5,  5,  5,  5, 10, 10, 10,  6,  0,  4,  6, 10,  6,  4,  4,  5,  5, 12, 14, 10,  0,  0,  0,  0, 12, 12,  5,  3,  3, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  7,  7,  7,  7,  0,  6,  6,  8,  8,  8,  6,  3,  2,  1,  4,  4,  7,  7,  8,  8,  0,  8,  0,  0,  8, 10,  0,  0,  0,  8,  0,  6,  0, 12,  6,  0,  7,  7,  7,  8,  5,  6,  6,  4, 10,  8,  8,  6, 10,  0, 10,  0, 10, 12, 12, 14, 14, 14,  6,  6, 10,  6,  7,  7,  7,  7, 10, 10, 10, 10,  0,  5,  7,  8,  8,  5,  5,  7,  7, 12, 14, 10,  0,  0,  0,  0, 14, 14,  8,  5,  4, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  3,  3,  5,  3,  0,  2,  4,  6,  6,  5,  6,  4,  4,  4,  1,  2,  5,  4,  5,  4,  0,  4,  0,  0,  7,  6,  0,  0,  0,  4,  0,  6,  0, 12,  6,  0,  6,  6,  5,  5,  6,  6,  5,  2, 10,  8,  6,  6,  5,  0,  7,  0,  8, 12, 12, 14, 12, 12,  3,  5,  7,  5,  6,  6,  5,  5, 10, 10, 10,  5,  0,  3,  5,  3,  6,  7,  7,  6,  5, 12, 14, 10,  0,  0,  0,  0, 14, 14, 10,  8,  7, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  5,  4,  5,  4,  0,  4,  2,  7,  6,  7,  6,  4,  4,  4,  2,  1,  5,  5,  5,  5,  0,  5,  0,  0,  7,  6,  0,  0,  0,  4,  0,  5,  0, 12,  5,  0,  5,  6,  5,  5,  3,  5,  6,  3, 10,  8,  6,  4, 10,  0,  8,  0,  8, 12, 12, 14, 12, 12,  4,  4,  8,  7,  5,  5,  5,  5, 10, 10, 10,  6,  0,  4,  7,  5,  6,  7,  8,  5,  5, 12, 14, 10,  0,  0,  0,  0, 14, 14, 10,  8,  7, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  5,  3,  5,  0,  7,  7,  7,  7,  8, 10,  8,  8,  7,  5,  5,  1,  2,  5,  5,  0,  5,  0,  0,  4,  5,  0,  0,  0,  4,  0, 10,  0,  8, 12,  0,  5,  6,  7,  7,  5,  8,  8,  4,  6,  4,  8,  7,  7,  0,  6,  0,  6,  8, 10, 12, 10,  8,  6,  7,  4, 10,  5,  5,  7,  7,  6,  6,  6,  8,  0, 10,  3,  2,  6, 12, 12,  5,  7,  7,  8,  5,  0,  0,  0,  0,  6,  6, 12, 12, 12, 10,  8,  8,  8,  8,  8,  8,  8 },
      {  7,  6,  4,  6,  0,  7,  7,  7,  7,  8, 10,  8,  8,  7,  4,  5,  2,  1,  5,  6,  0,  6,  0,  0,  3,  7,  0,  0,  0,  5,  0, 12,  0,  7, 12,  0,  5,  6,  8,  8,  3, 10, 10,  5,  3,  4,  8,  8,  6,  0,  4,  0,  4, 10, 12, 14,  6,  7,  7,  8,  5, 12,  5,  5,  8,  8,  3,  3,  3, 10,  0, 10,  2,  2,  7, 12, 12,  5,  8,  6,  8,  3,  0,  0,  0,  0,  5,  5, 12, 12, 12, 10,  8,  8,  8,  8,  8,  8,  8 },
      {  3,  3,  4,  3,  0,  5,  7,  7,  7,  8, 10,  6,  6,  8,  5,  5,  5,  5,  1,  3,  0,  3,  0,  0,  5,  6,  0,  0,  0,  3,  0, 10,  0, 10, 10,  0,  4,  6,  7,  7,  8,  8, 10,  4,  8,  7,  8,  8, 10,  0, 10,  0,  6, 10, 12, 14, 12, 10,  4,  6,  8, 12,  4,  4,  7,  7,  8,  8,  8,  8,  0, 10,  8,  5,  5, 10, 10,  4,  7, 10, 12,  8,  0,  0,  0,  0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  2,  2,  4,  2,  0,  4,  6,  7,  5,  6,  8,  6,  6,  8,  4,  5,  5,  6,  3,  1,  0,  2,  0,  0,  6,  7,  0,  0,  0,  3,  0, 10,  0, 10, 10,  0,  3,  6,  7,  6,  8,  7,  8,  4,  7,  6,  8,  7, 10,  0, 10,  0,  6, 10, 12, 14, 12, 10,  4,  4,  6, 10,  3,  3,  6,  6,  8,  8,  8,  8,  0, 10, 10,  4,  5, 10, 10,  3,  6, 10, 12,  8,  0,  0,  0,  0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  2,  2,  4,  2,  0,  4,  6,  7,  5,  6,  8,  6,  6,  8,  4,  5,  5,  6,  3,  2,  0,  1,  0,  0,  6,  7,  0,  0,  0,  3,  0, 10,  0, 10, 10,  0,  3,  6,  7,  6,  8,  7,  8,  4,  7,  6,  8,  7, 10,  0, 10,  0,  6, 10, 12, 14, 12, 10,  4,  4,  6, 10,  3,  3,  6,  6,  8,  8,  8,  8,  0, 10, 10,  4,  5, 10, 10,  3,  6, 10, 12,  8,  0,  0,  0,  0, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  6,  5,  4,  5,  0,  7,  7,  7,  7,  8,  8,  7,  7,  8,  7,  7,  4,  3,  5,  6,  0,  6,  0,  0,  1,  3,  0,  0,  0,  4,  0, 10,  0,  6, 12,  0,  4,  4,  6,  5,  4,  5,  7,  7,  4,  4,  8,  6,  8,  0,  7,  0,  5,  8, 10, 12, 10,  7,  7,  7,  5, 12,  4,  4,  6,  6,  4,  4,  4, 10,  0, 12,  5,  4,  3, 12, 12,  4,  6,  6,  8,  4,  0,  0,  0,  0,  6,  6, 12, 12, 12,  6,  6,  6,  6,  6,  6,  6,  6 },
      {  7,  8,  8,  8,  0,  6,  6,  5,  6,  6, 10,  8,  8, 10,  6,  6,  5,  7,  6,  7,  0,  7,  0,  0,  3,  1,  0,  0,  0,  8,  0, 10,  0,  8, 10,  0,  7,  7,  6,  7, 12,  4,  7,  7, 10,  8,  5, 10, 12,  0, 12,  0,  7,  6,  7,  7, 10,  6,  6,  6,  8,  7,  7,  7,  6,  6,  8,  8,  8,  7,  0, 10, 10,  7,  4,  6,  8,  7,  6,  8,  6,  8,  0,  0,  0,  0, 12, 12,  8,  8,  8, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  3,  2,  2,  2,  0,  4,  4,  7,  5,  6,  8,  6,  6,  8,  4,  4,  4,  5,  3,  3,  0,  3,  0,  0,  4,  8,  0,  0,  0,  1,  0, 10,  0,  8, 10,  0,  2,  4,  5,  5,  6,  4,  7,  4,  6,  5,  7,  7, 10,  0,  8,  0,  4, 10, 12, 14, 12, 10,  4,  3,  5, 10,  2,  2,  4,  4,  7,  7,  7,  8,  0, 10,  8,  3,  4, 10, 10,  2,  4, 10, 12,  7,  0,  0,  0,  0,  8,  8, 10, 10, 10,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      { 10, 10, 10, 10,  0,  5,  4, 12, 10,  6,  3,  4,  4,  6,  6,  5, 10, 12, 10, 10,  0, 10,  0,  0, 10, 10,  0,  0,  0, 10,  0,  1,  0, 20,  2,  0, 10, 12,  8, 10, 10,  4,  6,  8, 16, 14,  6,  6, 20,  0, 20,  0, 12, 20, 20, 20, 20, 20,  7,  8, 14,  2, 10, 10,  8,  8, 16, 16, 16,  5,  0,  2, 10,  8, 10,  2,  8, 10,  8, 20, 20, 16,  0,  0,  0,  0, 20, 20,  4,  4,  4, 20, 20, 20, 20, 20, 20, 20, 20 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      { 10,  8,  8,  8,  0, 12, 12,  4,  5,  8, 12, 10, 10, 12, 12, 12,  8,  7, 10, 10,  0, 10,  0,  0,  6,  8,  0,  0,  0,  8,  0, 20,  0,  1, 20,  0, 14,  4,  7,  5, 16, 12, 12, 16,  6,  7, 10, 12,  3,  0,  4,  0,  4,  2,  2,  2,  4,  2, 12, 14,  6, 20, 14, 14, 16, 16,  6,  6,  6, 12,  0, 20, 16, 10,  8, 20, 20, 14, 16,  4,  2,  6,  0,  0,  0,  0,  4,  4, 20, 20, 20,  4,  4,  4,  4,  4,  4,  4,  4 },
      { 10, 10, 12, 10,  0,  5,  4, 12, 10,  6,  3,  4,  4,  6,  6,  5, 12, 12, 10, 10,  0, 10,  0,  0, 12, 10,  0,  0,  0, 10,  0,  2,  0, 20,  1,  0, 10, 12,  8, 10, 10,  4,  6,  8, 16, 14,  6,  8, 20,  0, 20,  0, 12, 20, 20, 20, 20, 20,  7, 10, 14,  2, 10, 10,  8,  8, 16, 16, 16,  5,  0,  2, 10,  8, 10,  2,  8, 10,  8, 20, 20, 16,  0,  0,  0,  0, 20, 20,  4,  4,  4, 20, 20, 20, 20, 20, 20, 20, 20 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  3,  2,  2,  2,  0,  6,  6,  7,  5,  6,  7,  5,  5,  7,  6,  5,  5,  5,  4,  3,  0,  3,  0,  0,  4,  7,  0,  0,  0,  2,  0, 10,  0, 14, 10,  0,  1,  3,  5,  5,  5,  5,  7,  4,  7,  5,  6,  8, 14,  0,  8,  0,  7,  8, 10, 12, 10,  8,  4,  7,  7, 10,  1,  1,  5,  5,  7,  7,  7, 12,  0, 10,  6,  6,  7, 10, 10,  1,  5,  7, 12,  7,  0,  0,  0,  0, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  4,  4,  4,  0,  6,  6,  2,  2,  3,  7,  5,  5,  7,  6,  6,  6,  6,  6,  6,  0,  6,  0,  0,  4,  7,  0,  0,  0,  4,  0, 12,  0,  4, 12,  0,  3,  1,  5,  3,  8,  3,  8,  6,  7,  6,  7,  7,  7,  0,  6,  0,  4,  4,  5,  8,  8,  2,  6,  8,  6, 12,  3,  3,  6,  6,  6,  6,  6,  6,  0, 12,  8,  8,  6, 12, 12,  3,  5,  8,  6,  6,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  7,  5,  6,  5,  0,  4,  4,  4,  4,  3,  7,  5,  5,  7,  5,  5,  7,  8,  7,  7,  0,  7,  0,  0,  6,  6,  0,  0,  0,  5,  0,  8,  0,  7,  8,  0,  5,  5,  1,  3, 10,  5,  5,  4,  8,  7,  5, 10, 12,  0, 12,  0, 10, 12, 12, 14, 14, 12,  4,  8,  8,  8,  4,  4,  2,  2,  8,  8,  8, 10,  0, 10, 10, 10,  6,  8, 10, 12,  4,  7,  4,  8,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  6,  5,  6,  5,  0,  4,  4,  2,  2,  5,  8,  6,  6,  8,  5,  5,  7,  8,  7,  6,  0,  6,  0,  0,  5,  7,  0,  0,  0,  5,  0, 10,  0,  5, 10,  0,  5,  3,  3,  1, 10,  8,  8,  8, 10,  8,  5,  8, 10,  0, 10,  0,  8,  6,  5,  7,  8,  5,  8, 10,  8,  8,  7,  7,  5,  5,  7,  7,  7,  8,  0, 10, 10, 10,  8, 10, 12,  7,  5,  6,  5,  7,  0,  0,  0,  0, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  8,  6,  4,  6,  0,  6,  7, 12, 10,  8,  4,  5,  5,  5,  6,  3,  5,  3,  8,  8,  0,  8,  0,  0,  4, 12,  0,  0,  0,  6,  0, 10,  0, 16, 10,  0,  5,  8, 10, 10,  1,  8,  8, 10,  4,  5,  8, 10, 12,  0, 12,  0,  8, 12, 12, 14, 10, 12, 10,  8,  8, 10,  5,  5,  8,  8,  4,  4,  4,  8,  0,  8,  3,  6, 10,  4,  4,  5,  8, 10, 12,  4,  0,  0,  0,  0,  8,  8, 12, 12, 12,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  7,  4,  8,  4,  0,  4,  4,  3,  2,  3,  5,  4,  4,  6,  6,  5,  8, 10,  8,  7,  0,  7,  0,  0,  5,  4,  0,  0,  0,  4,  0,  4,  0, 12,  4,  0,  5,  3,  5,  8,  8,  1,  6,  7,  8,  7,  4,  6, 12,  0, 12,  0,  7,  8, 12, 12, 14, 12,  7,  5,  7,  8,  5,  5,  4,  4, 10, 10, 10,  5,  0,  5,  8,  8,  3,  4,  8,  5,  4, 14, 12, 10,  0,  0,  0,  0, 14, 14,  6,  6,  6, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  8,  7,  8,  7,  0,  4,  4,  6,  5,  2,  6,  5,  5,  6,  5,  6,  8, 10, 10,  8,  0,  8,  0,  0,  7,  7,  0,  0,  0,  7,  0,  6,  0, 12,  6,  0,  7,  8,  5,  8,  8,  6,  1,  8,  8,  8,  6,  8, 10,  0,  8,  0,  7,  7,  8, 10, 12,  7,  8,  8, 10,  3,  5,  5,  4,  4, 10, 10, 10,  6,  0,  7, 10, 10,  6,  4, 10,  5,  4, 12, 10, 10,  0,  0,  0,  0, 12, 12,  8,  8,  8, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  4,  4,  6,  4,  0,  4,  5,  7,  6,  6,  6,  4,  4,  4,  2,  3,  4,  5,  4,  4,  0,  4,  0,  0,  7,  7,  0,  0,  0,  4,  0,  8,  0, 16,  8,  0,  4,  6,  4,  8, 10,  7,  8,  1, 10,  8,  6,  8,  5,  0,  7,  0,  8, 10, 10, 12, 12, 10,  2,  7, 10,  8,  4,  4,  3,  3, 10, 10, 10,  6,  0,  6,  8,  8,  7,  7, 10,  4,  3, 12, 12, 10,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  7,  6,  6,  6,  0,  7, 10,  8,  8,  8, 10,  8,  8, 10, 10, 10,  6,  3,  8,  7,  0,  7,  0,  0,  4, 10,  0,  0,  0,  6,  0, 16,  0,  6, 16,  0,  7,  7,  8, 10,  4,  8,  8, 10,  1,  3, 10, 10,  4,  0,  4,  0,  4,  8, 10, 12, 12,  8, 10, 12,  8, 12,  7,  7,  8,  8,  4,  4,  4, 10,  0, 12,  4,  5,  7, 10,  8,  7,  8, 12, 12,  4,  0,  0,  0,  0,  8,  8, 12, 12, 12,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  6,  5,  3,  5,  0,  6,  8,  8,  7,  7,  8,  7,  7,  8,  8,  8,  4,  4,  7,  6,  0,  6,  0,  0,  4,  8,  0,  0,  0,  5,  0, 14,  0,  7, 14,  0,  5,  6,  7,  8,  5,  7,  8,  8,  3,  1,  8,  8,  8,  0,  7,  0,  6,  6,  8, 10, 10,  8,  6,  8,  6,  8,  5,  5,  8,  8,  3,  3,  3, 10,  0, 12,  3,  4,  7, 10,  8,  5,  8, 10, 10,  3,  0,  0,  0,  0,  8,  8, 14, 14, 14,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  8,  7,  8,  7,  0,  6,  6,  7,  5,  4,  8,  7,  7,  8,  6,  6,  8,  8,  8,  8,  0,  8,  0,  0,  8,  5,  0,  0,  0,  7,  0,  6,  0, 10,  6,  0,  6,  7,  5,  5,  8,  4,  6,  6, 10,  8,  1,  4, 10,  0, 10,  0,  5,  8,  8, 10, 12,  8,  6,  5,  8,  5,  6,  6,  5,  5,  8,  8,  8,  6,  0,  6,  8,  8,  5,  4,  8,  6,  5, 12, 10,  8,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  7,  6,  7,  6,  0,  6,  5,  7,  4,  5,  6,  6,  6,  6,  6,  4,  7,  8,  8,  7,  0,  7,  0,  0,  6, 10,  0,  0,  0,  7,  0,  6,  0, 12,  8,  0,  8,  7, 10,  8, 10,  6,  8,  8, 10,  8,  4,  1, 10,  0, 10,  0,  8, 10, 10, 12, 14, 10,  6,  7, 10,  7,  8,  8,  6,  6, 10, 10, 10,  8,  0,  6, 10, 10,  4,  7, 10,  8,  6, 14, 12, 10,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      { 10, 10,  8, 10,  0,  4, 10, 10, 10, 12, 12, 10, 10, 10,  5, 10,  7,  6, 10, 10,  0, 10,  0,  0,  8, 12,  0,  0,  0, 10,  0, 20,  0,  3, 20,  0, 14,  7, 12, 10, 12, 12, 10,  5,  4,  8, 10, 10,  1,  0,  2,  0,  8,  4,  4,  5,  6,  3,  4,  8,  6, 20, 14, 14, 16, 16,  4,  4,  4,  3,  0, 20,  8,  7, 10, 20, 16, 14, 16,  6,  5,  4,  0,  0,  0,  0,  6,  6, 20, 20, 20,  6,  6,  6,  6,  6,  6,  6,  6 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      { 10,  8,  6,  8,  0,  6, 10, 10,  8, 10, 12, 10, 10, 10,  7,  8,  6,  4, 10, 10,  0, 10,  0,  0,  7, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0,  8,  6, 12, 10, 12, 12,  8,  7,  4,  7, 10, 10,  2,  0,  1,  0,  6,  2,  5,  7,  8,  5,  8,  8, 12, 20,  8,  8, 10, 10,  7,  7,  7,  6,  0, 20, 10, 10,  8, 20, 16,  8, 10,  8,  7,  7,  0,  0,  0,  0,  6,  6, 20, 20, 20,  6,  6,  6,  6,  6,  6,  6,  6 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  6,  4,  6,  4,  0,  5,  5,  5,  4,  5, 10, 10, 10, 10,  8,  8,  6,  4,  6,  6,  0,  6,  0,  0,  5,  7,  0,  0,  0,  4,  0, 12,  0,  4, 12,  0,  7,  4, 10,  8,  8,  7,  7,  8,  4,  6,  5,  8,  8,  0,  6,  0,  1,  4,  6,  8, 10,  4, 10,  8,  5, 14,  7,  7,  8,  8,  4,  4,  4, 10,  0, 12,  6,  6,  8, 12, 10,  7,  8, 10,  8,  4,  0,  0,  0,  0,  8,  8, 14, 14, 14,  8,  8,  8,  8,  8,  8,  8,  8 },
      { 10, 10,  8, 10,  0, 12, 12,  4,  3,  6, 12, 12, 12, 12, 12, 12,  8, 10, 10, 10,  0, 10,  0,  0,  8,  6,  0,  0,  0, 10,  0, 20,  0,  2, 20,  0,  8,  4, 12,  6, 12,  8,  7, 10,  8,  6,  8, 10,  4,  0,  2,  0,  4,  1,  2,  2,  4,  2, 10, 10,  6, 20,  8,  8, 10, 10,  4,  4,  4, 10,  0, 20, 10, 10,  8, 20, 16,  8, 10,  4,  2,  4,  0,  0,  0,  0,  5,  5, 20, 20, 20,  5,  5,  5,  5,  5,  5,  5,  5 },
      { 12, 12, 10, 12,  0, 14, 14,  5,  3,  6, 12, 12, 12, 12, 12, 12, 10, 12, 12, 12,  0, 12,  0,  0, 10,  7,  0,  0,  0, 12,  0, 20,  0,  2, 20,  0, 10,  5, 12,  5, 12, 12,  8, 10, 10,  8,  8, 10,  4,  0,  5,  0,  6,  2,  1,  2,  4,  2, 10, 10,  6, 20,  8,  8, 10, 10,  4,  4,  4, 10,  0, 20, 10, 10,  8, 20, 16,  8, 10,  4,  2,  4,  0,  0,  0,  0,  5,  5, 20, 20, 20,  5,  5,  5,  5,  5,  5,  5,  5 },
      { 14, 14, 12, 14,  0, 16, 16,  5,  3,  8, 14, 14, 14, 14, 14, 14, 12, 14, 14, 14,  0, 14,  0,  0, 12,  7,  0,  0,  0, 14,  0, 20,  0,  2, 20,  0, 12,  8, 14,  7, 14, 12, 10, 12, 12, 10, 10, 12,  5,  0,  7,  0,  8,  2,  2,  1,  3,  3, 12, 12,  8, 20, 10, 10, 12, 12,  8,  8,  8, 12,  0, 20, 14, 14, 12, 20, 20, 10, 12,  4,  2,  8,  0,  0,  0,  0,  8,  8, 20, 20, 20,  8,  8,  8,  8,  8,  8,  8,  8 },
      { 12, 12, 10, 12,  0, 14, 14,  8,  6, 10, 12, 12, 12, 14, 12, 12, 10,  6, 12, 12,  0, 12,  0,  0, 10, 10,  0,  0,  0, 12,  0, 20,  0,  4, 20,  0, 10,  8, 14,  8, 10, 14, 12, 12, 12, 10, 12, 14,  6,  0,  8,  0, 10,  4,  4,  3,  1,  5, 12, 12, 10, 20,  8,  8, 14, 14, 10, 10, 10, 12,  0, 20, 12, 12, 12, 20, 20,  8, 14,  2,  4, 10,  0,  0,  0,  0,  8,  8, 20, 20, 20,  8,  8,  8,  8,  8,  8,  8,  8 },
      { 10, 10,  8, 10,  0, 12, 12,  4,  3,  5, 12, 12, 12, 14, 12, 12,  8,  7, 10, 10,  0, 10,  0,  0,  7,  6,  0,  0,  0, 10,  0, 20,  0,  2, 20,  0,  8,  2, 12,  5, 12, 12,  7, 10,  8,  8,  8, 10,  3,  0,  5,  0,  4,  2,  2,  3,  5,  1, 10, 12,  6, 20,  8,  8, 12, 12,  8,  8,  8, 12,  0, 20, 12, 12, 10, 20, 20,  8, 12,  5,  3,  8,  0,  0,  0,  0, 10, 10, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  4,  4,  6,  4,  0,  2,  3,  6,  5,  6,  5,  4,  4,  6,  3,  4,  6,  7,  4,  4,  0,  4,  0,  0,  7,  6,  0,  0,  0,  4,  0,  7,  0, 12,  7,  0,  4,  6,  4,  8, 10,  7,  8,  2, 10,  6,  6,  6,  4,  0,  8,  0, 10, 10, 10, 12, 12, 10,  1,  7, 10,  8,  4,  4,  3,  3, 10, 10, 10,  4,  0,  7, 10, 10,  7,  7, 10,  4,  3, 12, 12, 10,  0,  0,  0,  0, 14, 14,  7,  7,  7, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  4,  3,  6,  3,  0,  4,  4,  6,  5,  6,  6,  5,  5,  6,  5,  4,  7,  8,  6,  4,  0,  4,  0,  0,  7,  6,  0,  0,  0,  3,  0,  8,  0, 14, 10,  0,  7,  8,  8, 10,  8,  5,  8,  7, 12,  8,  5,  7,  8,  0,  8,  0,  8, 10, 10, 12, 12, 12,  7,  1,  5,  7,  7,  7,  7,  7, 12, 12, 12,  8,  0,  8, 10, 10,  7,  6, 10,  7,  7,  8, 12, 12,  0,  0,  0,  0, 14, 14,  8,  8,  8, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  5,  3,  5,  0,  8,  8,  6,  5,  6,  8,  8,  8, 10,  7,  8,  4,  5,  8,  6,  0,  6,  0,  0,  5,  8,  0,  0,  0,  5,  0, 14,  0,  6, 14,  0,  7,  6,  8,  8,  8,  7, 10, 10,  8,  6,  8, 10,  6,  0, 12,  0,  5,  6,  6,  8, 10,  6, 10,  5,  1, 10,  7,  7,  8,  8,  3,  3,  3,  8,  0, 14,  5,  5, 10, 10, 10,  7,  8, 10,  8, 10,  0,  0,  0,  0,  8,  8, 16, 16, 16,  8,  8,  8,  8,  8,  8,  8,  8 },
      { 10, 10, 12, 10,  0,  6,  6,  5,  4,  3,  5,  6,  6,  6,  5,  7, 10, 12, 12, 10,  0, 10,  0,  0, 12,  7,  0,  0,  0, 10,  0,  2,  0, 20,  2,  0, 10, 12,  8,  8, 10,  8,  3,  8, 12,  8,  5,  7, 20,  0, 20,  0, 14, 20, 20, 20, 20, 20,  8,  7, 10,  1, 10, 10,  8,  8, 20, 20, 20,  5,  0,  2, 10, 10,  7,  2,  6, 10,  8, 20, 20, 20,  0,  0,  0,  0, 20, 20,  4,  4,  4, 20, 20, 20, 20, 20, 20, 20, 20 },
      {  3,  2,  2,  2,  0,  6,  6,  7,  5,  6,  7,  5,  5,  7,  6,  5,  5,  5,  4,  3,  0,  3,  0,  0,  4,  7,  0,  0,  0,  2,  0, 10,  0, 14, 10,  0,  1,  3,  4,  7,  5,  5,  5,  4,  7,  5,  6,  8, 14,  0,  8,  0,  7,  8,  8, 10,  8,  8,  4,  7,  7, 10,  1,  2,  5,  5,  7,  7,  7, 12,  0, 10,  6,  6,  7, 10, 10,  2,  5, 10, 12,  7,  0,  0,  0,  0, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  3,  2,  2,  2,  0,  6,  6,  7,  5,  6,  7,  5,  5,  7,  6,  5,  5,  5,  4,  3,  0,  3,  0,  0,  4,  7,  0,  0,  0,  2,  0, 10,  0, 14, 10,  0,  1,  3,  4,  7,  5,  5,  5,  4,  7,  5,  6,  8, 14,  0,  8,  0,  7,  8,  8, 10,  8,  8,  4,  7,  7, 10,  2,  1,  5,  5,  7,  7,  7, 12,  0, 10,  6,  6,  7, 10, 10,  2,  5, 10, 12,  7,  0,  0,  0,  0, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  4,  6,  4,  0,  4,  4,  4,  4,  4,  6,  5,  5,  7,  5,  5,  7,  8,  7,  6,  0,  6,  0,  0,  6,  6,  0,  0,  0,  4,  0,  8,  0, 16,  8,  0,  5,  6,  2,  5,  8,  4,  4,  3,  8,  8,  5,  6, 16,  0, 10,  0,  8, 10, 10, 12, 14, 12,  3,  7,  8,  8,  5,  5,  1,  2,  6,  6,  6, 10,  0,  8,  6,  6,  6,  8, 10,  5,  2,  6,  4,  6,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  6,  4,  6,  4,  0,  4,  4,  4,  4,  4,  6,  5,  5,  7,  5,  5,  7,  8,  7,  6,  0,  6,  0,  0,  6,  6,  0,  0,  0,  4,  0,  8,  0, 16,  8,  0,  5,  6,  2,  5,  8,  4,  4,  3,  8,  8,  5,  6, 16,  0, 10,  0,  8, 10, 10, 12, 14, 12,  3,  7,  8,  8,  5,  5,  2,  1,  6,  6,  6, 10,  0,  8,  6,  6,  6,  8, 10,  5,  2,  6,  4,  6,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      {  8,  7,  4,  7,  0,  8,  8,  7,  7,  7, 10, 10, 10, 10, 10, 10,  6,  3,  8,  8,  0,  8,  0,  0,  4,  8,  0,  0,  0,  7,  0, 16,  0,  6, 16,  0,  7,  6,  8,  7,  4, 10, 10, 10,  4,  3,  8, 10,  4,  0,  7,  0,  4,  4,  4,  8, 10,  8, 10, 12,  3, 20,  7,  7,  6,  6,  1,  2,  2, 10,  0, 20,  5,  5, 10, 20, 16,  7,  6,  8,  6,  2,  0,  0,  0,  0, 10, 10, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  8,  7,  4,  7,  0,  8,  8,  7,  7,  7, 10, 10, 10, 10, 10, 10,  6,  3,  8,  8,  0,  8,  0,  0,  4,  8,  0,  0,  0,  7,  0, 16,  0,  6, 16,  0,  7,  6,  8,  7,  4, 10, 10, 10,  4,  3,  8, 10,  4,  0,  7,  0,  4,  4,  4,  8, 10,  8, 10, 12,  3, 20,  7,  7,  6,  6,  2,  1,  2, 10,  0, 20,  5,  5, 10, 20, 16,  7,  6,  8,  6,  2,  0,  0,  0,  0, 10, 10, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  8,  7,  4,  7,  0,  8,  8,  7,  7,  7, 10, 10, 10, 10, 10, 10,  6,  3,  8,  8,  0,  8,  0,  0,  4,  8,  0,  0,  0,  7,  0, 16,  0,  6, 16,  0,  7,  6,  8,  7,  4, 10, 10, 10,  4,  3,  8, 10,  4,  0,  7,  0,  4,  4,  4,  8, 10,  8, 10, 12,  3, 20,  7,  7,  6,  6,  2,  2,  1, 10,  0, 20,  5,  5, 10, 20, 16,  7,  6,  8,  6,  2,  0,  0,  0,  0, 10, 10, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  8,  8, 10,  8,  0,  3,  4,  6,  7,  8,  5,  6,  6, 10,  5,  6,  8, 10,  8,  8,  0,  8,  0,  0, 10,  7,  0,  0,  0,  8,  0,  5,  0, 12,  5,  0, 12,  6, 10,  8,  8,  5,  6,  6, 10, 10,  6,  8,  3,  0,  6,  0, 10, 10, 10, 12, 12, 12,  4,  8,  8,  5, 12, 12, 10, 10, 10, 10, 10,  1,  0,  4,  7,  7,  6,  5,  8, 12, 10, 12, 14, 10,  0,  0,  0,  0, 10, 10,  5,  4,  4, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      { 10, 10, 12, 10,  0,  6,  6, 10, 10,  8,  4,  4,  4,  5,  3,  4, 10, 10, 10, 10,  0, 10,  0,  0, 12, 10,  0,  0,  0, 10,  0,  2,  0, 20,  2,  0, 10, 12, 10, 10,  8,  5,  7,  6, 12, 12,  6,  6, 20,  0, 20,  0, 12, 20, 20, 20, 20, 20,  7,  8, 14,  2, 10, 10,  8,  8, 20, 20, 20,  4,  0,  1,  6,  6,  8,  2,  6, 10, 10, 20, 20, 20,  0,  0,  0,  0, 20, 20,  4,  3,  3, 20, 20, 20, 20, 20, 20, 20, 20 },
      { 10,  8,  5,  8,  0,  6,  6, 10, 10,  8,  6,  6,  6,  7,  5,  7,  3,  2,  8, 10,  0, 10,  0,  0,  5, 10,  0,  0,  0,  8,  0, 10,  0, 16, 10,  0,  6,  8, 10, 10,  3,  8, 10,  8,  4,  3,  8, 10,  8,  0, 10,  0,  6, 10, 10, 14, 12, 12, 10, 10,  5, 10,  6,  6,  6,  6,  5,  5,  5,  7,  0,  6,  1,  3,  8,  7,  5,  6, 10, 12, 14,  5,  0,  0,  0,  0,  8,  8, 14, 14, 14,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  4,  3,  3,  3,  0,  7,  8, 10, 10, 10, 10, 10, 10,  8,  3,  5,  2,  2,  5,  4,  0,  4,  0,  0,  4,  7,  0,  0,  0,  3,  0,  8,  0, 10,  8,  0,  6,  8, 10, 10,  6,  8, 10,  8,  5,  4,  8, 10,  7,  0, 10,  0,  6, 10, 10, 14, 12, 12, 10, 10,  5, 10,  6,  6,  6,  6,  5,  5,  5,  7,  0,  6,  3,  1,  7, 10,  7,  6, 10,  7,  8,  5,  0,  0,  0,  0,  8,  8, 14, 14, 14,  8,  8,  8,  8,  8,  8,  8,  8 },
      {  5,  4,  6,  4,  0,  6,  5,  6,  3,  6,  6,  6,  6,  8,  6,  6,  6,  7,  5,  5,  0,  5,  0,  0,  3,  4,  0,  0,  0,  4,  0, 10,  0,  8, 10,  0,  7,  6,  6,  8, 10,  3,  6,  7,  7,  7,  5,  4, 10,  0,  8,  0,  8,  8,  8, 12, 12, 10,  7,  7, 10,  7,  7,  7,  6,  6, 10, 10, 10,  6,  0,  8,  8,  7,  1,  7,  8,  7,  6, 10,  8, 10,  0,  0,  0,  0, 10, 10,  7,  7,  7, 10, 10, 10, 10, 10, 10, 10, 10 },
      { 10, 10, 12, 10,  0,  6,  6,  6,  7,  4,  3,  4,  4,  5,  7,  7, 12, 12, 10, 10,  0, 10,  0,  0, 12,  6,  0,  0,  0, 10,  0,  2,  0, 20,  2,  0, 10, 12,  8, 10,  4,  4,  4,  7, 10, 10,  4,  7, 20,  0, 20,  0, 12, 20, 20, 20, 20, 20,  7,  6, 10,  2, 10, 10,  8,  8, 20, 20, 20,  5,  0,  2,  7, 10,  7,  1,  8, 10,  8, 20, 20, 20,  0,  0,  0,  0, 20, 20,  4,  4,  4, 20, 20, 20, 20, 20, 20, 20, 20 },
      { 10, 10, 12, 10,  0,  8,  8, 10, 10, 10,  5,  4,  4,  5,  7,  8, 12, 12, 10, 10,  0, 10,  0,  0, 12,  8,  0,  0,  0, 10,  0,  8,  0, 20,  8,  0, 10, 12, 10, 12,  4,  8, 10, 10,  8,  8,  8, 10, 16,  0, 16,  0, 10, 16, 16, 20, 20, 20, 10, 10, 10,  6, 10, 10, 10, 10, 16, 16, 16,  8,  0,  6,  5,  7,  8,  8,  1, 10, 10, 20, 20, 16,  0,  0,  0,  0, 10, 10,  8,  6,  7, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  3,  2,  2,  2,  0,  5,  6,  7,  5,  6,  7,  5,  5,  7,  6,  5,  5,  5,  4,  3,  0,  3,  0,  0,  4,  7,  0,  0,  0,  2,  0, 10,  0, 14, 10,  0,  1,  3, 12,  7,  5,  5,  5,  4,  7,  5,  6,  8, 14,  0,  8,  0,  7,  8,  8, 10,  8,  8,  4,  7,  7, 10,  2,  2,  5,  5,  7,  7,  7, 12,  0, 10,  6,  6,  7, 10, 10,  1,  5,  3,  6,  7,  0,  0,  0,  0, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14, 14 },
      {  6,  4,  4,  4,  0,  4,  6,  4,  4,  4,  6,  5,  5,  7,  5,  5,  7,  8,  7,  6,  0,  6,  0,  0,  6,  6,  0,  0,  0,  4,  0,  8,  0, 16,  8,  0,  5,  5,  4,  5,  8,  4,  4,  3,  8,  8,  5,  6, 16,  0, 10,  0,  8, 10, 10, 12, 14, 12,  3,  7,  8,  8,  5,  5,  2,  2,  6,  6,  6, 10,  0, 10, 10, 10,  6,  8, 10,  5,  1,  6,  4,  6,  0,  0,  0,  0, 12, 12, 10, 10, 10, 12, 12, 12, 12, 12, 12, 12, 12 },
      { 10, 10,  7, 10,  0, 12, 12,  7,  6,  8, 12, 12, 12, 12, 12, 12,  7,  6, 10, 10,  0, 10,  0,  0,  6,  8,  0,  0,  0, 10,  0, 20,  0,  4, 20,  0,  7,  8,  7,  6, 10, 14, 12, 12, 12, 10, 12, 14,  6,  0,  8,  0, 10,  4,  4,  4,  2,  5, 12,  8, 10, 20, 10, 10,  6,  6,  8,  8,  8, 12,  0, 20, 12,  7, 10, 20, 20,  3,  6,  1,  3,  8,  0,  0,  0,  0,  8,  8, 20, 20, 20,  8,  8,  8,  8,  8,  8,  8,  8 },
      { 12, 12, 10, 12,  0, 10, 10,  5,  4,  6, 14, 14, 14, 14, 14, 14,  8,  8, 12, 12,  0, 12,  0,  0,  8,  6,  0,  0,  0, 12,  0, 20,  0,  2, 20,  0, 12,  6,  4,  5, 12, 12, 10, 12, 12, 10, 10, 12,  5,  0,  7,  0,  8,  2,  2,  2,  4,  3, 12, 12,  8, 20, 12, 12,  4,  4,  6,  6,  6, 14,  0, 20, 14,  8,  8, 20, 20,  6,  4,  3,  1,  6,  0,  0,  0,  0,  5,  5, 20, 20, 20,  5,  5,  5,  5,  5,  5,  5,  5 },
      {  8,  7,  5,  7,  0,  8,  8,  7,  7,  7, 10, 10, 10, 10, 10, 10,  5,  3,  8,  8,  0,  8,  0,  0,  4,  8,  0,  0,  0,  7,  0, 16,  0,  6, 16,  0,  7,  6,  8,  7,  4, 10, 10, 10,  4,  3,  8, 10,  4,  0,  7,  0,  4,  4,  4,  8, 10,  8, 10, 12, 10, 20,  7,  7,  6,  6,  2,  2,  2, 10,  0, 20,  5,  5, 10, 20, 16,  7,  6,  8,  6,  1,  0,  0,  0,  0, 10, 10, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      {  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  1,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0,  0 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  6,  5, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  1,  2, 20, 20, 20,  6,  5,  4,  4,  4,  4,  4,  4 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  6,  5, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  2,  1, 20, 20, 20,  7,  6,  5,  5,  5,  5,  5,  5 },
      { 10, 10, 12, 10,  0,  8,  8, 10, 10, 10,  3,  4,  5,  8, 10, 10, 12, 12, 10, 10,  0, 10,  0,  0, 12,  8,  0,  0,  0, 10,  0,  4,  0, 20,  4,  0, 14, 10, 10, 12, 12,  6,  8,  8, 12, 14,  8,  8, 20,  0, 20,  0, 14, 20, 20, 20, 20, 20,  7,  8, 16,  4, 14, 14, 10, 10, 20, 20, 20,  5,  0,  4, 14, 14,  7,  4,  8, 14, 10, 20, 20, 20,  0,  0,  0,  0, 20, 20,  1,  3,  4, 20, 20, 20, 20, 20, 20, 20, 20 },
      { 10, 10, 12, 10,  0,  8,  8, 10, 10, 10,  4,  3,  3,  5,  8,  8, 12, 12, 10, 10,  0, 10,  0,  0, 12,  8,  0,  0,  0, 10,  0,  4,  0, 20,  4,  0, 14, 10, 10, 12, 12,  6,  8,  8, 12, 14,  8,  8, 20,  0, 20,  0, 14, 20, 20, 20, 20, 20,  7,  8, 16,  4, 14, 14, 10, 10, 20, 20, 20,  4,  0,  3, 14, 14,  7,  4,  6, 14, 10, 20, 20, 20,  0,  0,  0,  0, 20, 20,  3,  1,  2, 20, 20, 20, 20, 20, 20, 20, 20 },
      { 10, 10, 12, 10,  0,  8,  8, 10, 10, 10,  5,  3,  3,  4,  7,  7, 12, 12, 10, 10,  0, 10,  0,  0, 12,  8,  0,  0,  0, 10,  0,  4,  0, 20,  4,  0, 14, 10, 10, 12, 12,  6,  8,  8, 12, 14,  8,  8, 20,  0, 20,  0, 14, 20, 20, 20, 20, 20,  7,  8, 16,  4, 14, 14, 10, 10, 20, 20, 20,  4,  0,  3, 14, 14,  7,  4,  7, 14, 10, 20, 20, 20,  0,  0,  0,  0, 20, 20,  4,  2,  1, 20, 20, 20, 20, 20, 20, 20, 20 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14, 10, 10, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  6,  7, 20, 20, 20,  1,  2,  3,  5,  5,  5,  5,  5 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  5,  6, 20, 20, 20,  2,  1,  2,  4,  4,  4,  4,  4 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  3,  2,  1,  4,  4,  4,  4,  4 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  5,  4,  4,  1,  3,  3,  3,  3 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  5,  4,  4,  3,  1,  5,  4,  2 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  5,  4,  4,  3,  5,  1,  3,  4 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  5,  4,  4,  3,  4,  3,  1,  3 },
      { 10,  8,  6,  8,  0, 14, 14, 12, 12, 12, 14, 12, 12, 14, 14, 14,  8,  8, 10, 10,  0, 10,  0,  0,  6, 12,  0,  0,  0,  8,  0, 20,  0,  4, 20,  0, 14, 12, 12, 12,  8, 14, 12, 14,  8,  8, 14, 14,  6,  0,  6,  0,  8,  5,  5,  8,  8, 10, 14, 14,  8, 20, 14, 14, 12, 12, 10, 10, 10, 10,  0, 20,  8,  8, 10, 20, 10, 14, 12,  8,  5, 10,  0,  0,  0,  0,  4,  5, 20, 20, 20,  5,  4,  4,  3,  2,  4,  3,  1 },
    };


    MolHistogram::MolHistogram(const ROMol &mol, const double *dmat) :
      d_h(boost::extents[mol.getNumHeavyAtoms()][O3_MAX_H_BINS]) {
      PRECONDITION(dmat,"empty distance matrix");
      unsigned int nAtoms = mol.getNumAtoms();
      unsigned int y = 0;
      for (unsigned int i = 0; i < nAtoms; ++i) {
        // skip hydrogens
        if ((mol.getAtomWithIdx(i))->getAtomicNum() == 1) {
          continue;
        }
        // initialize h row(y) to all zeros
        for (unsigned int j = 0; j < O3_MAX_H_BINS; ++j) {
          d_h[y][j] = 0;
        }
        for (unsigned int j = 0; j < nAtoms; ++j) {
          unsigned int dist = static_cast<unsigned int>(floor(dmat[i * nAtoms + j]));
          if (dist < O3_MAX_H_BINS) {
            ++d_h[y][dist];
          }
        }
        ++y;
      }
    }


    int o3aMMFFCostFunc(const unsigned int prbIdx,
      const unsigned int refIdx, double hSum, void *data)
    {
      boost::uint8_t mmffSim = mmffSimMatrix
        [(static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFAtomType(refIdx) - 1]
        [(static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFAtomType(prbIdx) - 1];
      
      return boost::math::iround
        ((static_cast<double>(((O3AFuncData *)data)->coeff) * O3_CHARGE_WEIGHT
        * fabs((static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFPartialCharge(refIdx)
        - (static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFPartialCharge(prbIdx))
        + (((O3AFuncData *)data)->useMMFFSim
        ? static_cast<double>(O3_MAX_WEIGHT_COEFF - ((O3AFuncData *)data)->coeff)
        * static_cast<double>(mmffSim) : 0.0) + hSum) * 1.0e03);
    }


    int o3aCrippenCostFunc(const unsigned int prbIdx,
      const unsigned int refIdx, double hSum, void *data)
    {
      return boost::math::iround
        ((static_cast<double>((static_cast<O3AFuncData *>(data))->coeff)
        * O3_CHARGE_WEIGHT * fabs((*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->refProp)))[refIdx]
        - (*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->prbProp)))[prbIdx])
        + hSum) * 1.0e03);
    }


    void LAP::computeCostMatrix(const ROMol &prbMol, const MolHistogram &prbHist,
      const ROMol &refMol, const MolHistogram &refHist, O3AConstraintVect *o3aConstraintVect,
      int (*costFunc)(const unsigned int, const unsigned int, double, void *),
      void *data, const unsigned int n_bins)
    {
      unsigned int i;
      unsigned int j;
      unsigned int k;
      unsigned int x;
      unsigned int y;
      unsigned int largestNHeavyAtoms =
        std::max(refMol.getNumHeavyAtoms(), prbMol.getNumHeavyAtoms());
      unsigned int mapIdx = 0;
      double hSum = 0.0;
      
      
      // ref is on rows, prb is on columns
      for (i = 0; i < largestNHeavyAtoms; ++i) {
        for (j = 0; j < largestNHeavyAtoms; ++j) {
          d_cost[i][j] = O3_DUMMY_COST;
        }
      }
      for (i = 0, y = 0; i < refMol.getNumAtoms(); ++i) {
        if (refMol[i]->getAtomicNum() == 1) {
          continue;
        }
        for (j = 0, x = 0; j < prbMol.getNumAtoms(); ++j) {
          if (prbMol[j]->getAtomicNum() == 1) {
            continue;
          }
          // if this pair is constrained, make sure that it is
          // associated with a low cost so it is picked up
          // in the minimum cost path
          // first pair element is prb, second is ref
          if (o3aConstraintVect && (mapIdx < o3aConstraintVect->size())
            && (j == (*o3aConstraintVect)[mapIdx]->getPrbIdx())
            && (i == (*o3aConstraintVect)[mapIdx]->getRefIdx())) {
            d_cost[y][x] = O3_LARGE_NEGATIVE_WEIGHT;
            ++mapIdx;
          }
          else {
            for (k = 0, hSum = 0.0; k < n_bins; ++k) {
              int rhyk = refHist.get(y,k);
              int phxk = prbHist.get(x,k);
              if ((!rhyk) && (!phxk)) continue;
              hSum += (double)square(rhyk - phxk) / (double)(rhyk + phxk);
            }
            d_cost[y][x] = (*costFunc)(j, i, hSum, data);
          }
          ++x;
        }
        ++y;
      }
    }


    //----------------------------------------------------
    // Taken from Jonker, R.; Volgenant, A.
    // A Shortest Augmenting Path Algorithm for Dense
    // and Sparse Linear Assignment Problems
    // Computing 1987, 38, 325-340
    //----------------------------------------------------
    void LAP::computeMinCostPath(const int dim)
    {
      // input:
      // dim    - problem size
      // cost   - cost matrix

      // output:
      // rowsol - column assigned to row in solution
      // colsol - row assigned to column in solution
      // v      - dual variables, column reduction numbers

      int loopCnt;
      int unassignedFound;
      // row vars
      int i;
      int imin;
      int numFree = 0;
      int prvNumFree;
      int f;
      int i0;
      int k;
      int freeRow;
      // col vars
      int j;
      int j1;
      int j2 = 0;
      int endOfPath;
      int last = 0;
      int low;
      int up;
      // cost vars
      int min = 0;
      int h;
      int uMin;
      int uSubMin;
      int v2;

      // init how many times a row will be assigned in the column reduction
      for (int n = 0; n < dim; ++n) {
        d_matches[n] = 0;
      }
      // COLUMN REDUCTION
      // reverse order gives better results
      for (j = dim - 1; j >= 0; --j) {
        // find minimum cost over rows
        min = d_cost[0][j];
        imin = 0;
        for (i = 1; i < dim; ++i) {
          if (d_cost[i][j] < min) {
            min = d_cost[i][j];
            imin = i;
          }
        }
        d_v[j] = min;

        if (++(d_matches[imin]) == 1) {
          // init assignment if minimum row assigned for first time
          d_rowSol[imin] = j;
          d_colSol[j] = imin;
        }
        else {
          // row already assigned, column not assigned
          d_colSol[j] = -1;
        }
      }
      // REDUCTION TRANSFER
      for (i = 0; i < dim; ++i) {
        if (!(d_matches[i])) {
          // fill list of unassigned 'd_free' rows
          d_free[numFree++] = i;
        }
        else {
          if (d_matches[i] == 1) {
            // transfer reduction from rows that are assigned once
            j1 = d_rowSol[i];
            min = O3_DUMMY_COST;
            for (j = 0; j < dim; ++j) {
              if ((j != j1) && (d_cost[i][j] - d_v[j] < min)) {
                min = d_cost[i][j] - d_v[j];
              }
            }
            d_v[j1] -= min;
          }
        }
      }
      // AUGMENTING ROW REDUCTION
      // do-loop to be done twice
      for (loopCnt = 0; loopCnt < 2; ++loopCnt) {
        // scan all d_free rows
        // in some cases, a d_free row may be replaced
        // with another one to be scanned next
        k = 0;
        prvNumFree = numFree;
        numFree = 0;
        // start list of rows still d_free after augmenting row reduction
        while (k < prvNumFree) {
          i = d_free[k];
          ++k;
          // find minimum and second minimum reduced cost over columns
          uMin = d_cost[i][0] - d_v[0];
          j1 = 0;
          uSubMin = O3_DUMMY_COST;
          for (j = 1; j < dim; ++j) {
            h = d_cost[i][j] - d_v[j];
            if (h < uSubMin) {
              if (h >= uMin) {
                uSubMin = h;
                j2 = j;
              }
              else {
                uSubMin = uMin;
                uMin = h;
                j2 = j1;
                j1 = j;
              }
            }
          }
          i0 = d_colSol[j1];
          if (uMin < uSubMin) {
            // change the reduction of the minimum column to increase
            // the minimum reduced cost in the row to the subminimum
            d_v[j1] -= (uSubMin - uMin);
          }
          else {
            // minimum and subminimum equal
            if (i0 >= 0) {
              // minimum column j1 is assigned
              // swap columns j1 and j2, as j2 may be unassigned
              j1 = j2;
              i0 = d_colSol[j2];
            }
          }
          // (re-)assign i to j1, possibly de-assigning an i0
          d_rowSol[i] = j1;
          d_colSol[j1] = i;

          if (i0 >= 0) {
            // minimum column j1 assigned earlier
            if (uMin < uSubMin) {
              // put in current k, and go back to that k
              // continue augmenting path i - j1 with i0
              d_free[--k] = i0;
            }
            else {
              // no further augmenting reduction possible
              // store i0 in list of d_free rows for next phase
              d_free[numFree++] = i0;
            }
          }
        }
      }
      // AUGMENT SOLUTION for each d_free row
      for (f = 0; f < numFree; ++f) {
        // start row of augmenting path
        freeRow = d_free[f];
        // Dijkstra shortest path algorithm
        // runs until unassigned column added to shortest path tree
        for (j = 0; j < dim; ++j) {
          d_d[j] = d_cost[freeRow][j] - d_v[j];
          d_pred[j] = freeRow;
          // init column list
          d_colList[j] = j;
        }
        // columns in 0..low-1 are ready, now none
        low = 0;
        // columns in low..up-1 are to be scanned for current minimum, now none
        up = 0;
        // columns in up..dim-1 are to be considered later to find new minimum,
        // at this stage the list simply contains all columns
        unassignedFound = 0;
        while (!unassignedFound) {
          if (up == low) {
            // no more columns to be scanned for current minimum
            last = low - 1;
            // scan columns for up..dim-1 to find all indices for which new minimum occurs
            // store these indices between low..up-1 (increasing up)
            min = d_d[d_colList[up++]];
            for (k = up; k < dim; ++k) {
              j = d_colList[k];
              h = d_d[j];
              if (h <= min) {
                if (h < min) {
                  // new minimum
                  // restart list at index low
                  up = low;
                  min = h;
                }
                // new index with same minimum, put on undex up, and extend list
                d_colList[k] = d_colList[up];
                d_colList[up++] = j;
              }
            }
            // check if any of the minimum columns happens to be unassigned
            // if so, we have an augmenting path right away
            for (k = low; k < up; ++k) {
              if (d_colSol[d_colList[k]] < 0) {
                endOfPath = d_colList[k];
                unassignedFound = 1;
                break;
              }
            }
          }
          if (!unassignedFound) {
            // update 'distances' between freeRow and
            // all unscanned columns, via next scanned column
            j1 = d_colList[low];
            ++low;
            i = d_colSol[j1];
            h = d_cost[i][j1] - d_v[j1] - min;
            for (k = up; k < dim; ++k) {
              j = d_colList[k];
              v2 = d_cost[i][j] - d_v[j] - h;
              if (v2 < d_d[j]) {
                d_pred[j] = i;
                if (v2 == min) {
                  // new column found at same minimum value
                  if (d_colSol[j] < 0) {
                    // if unassigned, shortest augmenting path is complete
                    endOfPath = j;
                    unassignedFound = 1;
                    break;
                  }
                  else {
                    // else add to list to be scanned right away
                    d_colList[k] = d_colList[up];
                    d_colList[up++] = j;
                  }
                }
                d_d[j] = v2;
              }
            }
          }
        }
        // update column prices
        for (k = 0; k <= last; ++k) {
          j1 = d_colList[k];
          d_v[j1] += (d_d[j1] - min);
        }
        // reset row and column assignments along the alternating path
        do {
          i = d_pred[endOfPath];
          d_colSol[endOfPath] = i;
          j1 = endOfPath;
          endOfPath = d_rowSol[i];
          d_rowSol[i] = j1;
        } while (i != freeRow);
      }
    }


    void SDM::fillFromLAP(LAP &lap)
    {
      unsigned int i;
      unsigned int j;
      unsigned int k;
      unsigned int n;
      unsigned int n_atoms;
      unsigned int n_equiv;
      const RDGeom::POINT3D_VECT &refPos = d_refConf->getPositions();
      const RDGeom::POINT3D_VECT &prbPos = d_prbConf->getPositions();
      const ROMol *mol[2] = { &(d_refConf->getOwningMol()), &(d_prbConf->getOwningMol()) };
      
      
      // filter out atom equivalences whose cost is > O3_DUMMY_COST
      for (i = 0, n_equiv = 0; i < mol[0]->getNumHeavyAtoms(); ++i) {
        if (lap.getCost(i, lap.getRowSol(i)) >= O3_DUMMY_COST) {
          continue;
        }
        SDMElement *sdmElement = new SDMElement;
        sdmElement->idx[0] = i;
        sdmElement->idx[1] = lap.getRowSol(i);
        sdmElement->cost = lap.getCost(i, lap.getRowSol(i));
        d_SDMPtrVect.push_back(boost::shared_ptr<SDMElement>(sdmElement));
        ++n_equiv;
      }
      boost::multi_array<double, 2> diff(boost::extents[n_equiv][n_equiv]);
      // find out correspondences between heavy atom indexes
      // in the original molecule and heavy atoms in the sdm matrix
      // (the sdm matrix is rebuilt using original molecule
      // heavy atom indexes)
      for (n = 0; n < 2; ++n) {
        for (i = 0; i < n_equiv; ++i) {
          j = 0;
          n_atoms = 0;
          k = d_SDMPtrVect[i]->idx[n];
          while ((n_atoms < mol[n]->getNumAtoms()) && (j <= k)) {
            if (mol[n]->getAtomWithIdx(n_atoms)->getAtomicNum() > 1) {
              ++j;
              d_SDMPtrVect[i]->idx[n] = n_atoms;
            }
            ++n_atoms;
          }
        }
      }
      // loop over n_equiv rows
      for (i = 0; i < n_equiv; ++i) {
        d_SDMPtrVect[i]->o3aConstraint = NULL;
        if (d_SDMPtrVect[i]->cost < 0) {
          // this is one of the constrained atom pairs, so assign it
          // a pointer to the relevant O3AConstraint object
          if (d_o3aConstraintVect) {
            for (j = 0; j < (*d_o3aConstraintVect).size(); ++j) {
              if (((*d_o3aConstraintVect)[j]->getPrbIdx() == d_SDMPtrVect[i]->idx[1])
                && ((*d_o3aConstraintVect)[j]->getRefIdx() == d_SDMPtrVect[i]->idx[0])) {
                d_SDMPtrVect[i]->o3aConstraint = (*d_o3aConstraintVect)[j];
                break;
              }
            }
          }
        }
        // initialize to 0 the i-th row of the diff matrix
        for (j = 0; j < n_equiv; ++j) {
          diff[i][j] = 0.0;
        }
        for (j = 0; j < n_equiv; ++j) {
          // if i == j then the distance will be 0,
          // no need to spend time computing it
          if (i == j) {
            continue;
          }
          // - the euclidean distance between atoms i,j is computed for ref_conf
          // - the euclidean distance between the corresponding atoms is computed for moved_conf
          // - the absolute value of the difference between the two distances is computed
          //   and placed in diff[i][j]
          diff[i][j] = fabs((refPos[d_SDMPtrVect[i]->idx[0]] - refPos[d_SDMPtrVect[j]->idx[0]]).length()
            - (prbPos[d_SDMPtrVect[i]->idx[1]] - prbPos[d_SDMPtrVect[j]->idx[1]]).length());
        }
      }
      for (i = 0; i < n_equiv; ++i) {
        for (j = 0, d_SDMPtrVect[i]->score = 0; j < n_equiv; ++j) {
          if ((diff[i][j] > O3_THRESHOLD_DIFF_DISTANCE) && (d_SDMPtrVect[i]->cost >= 0)) {
            ++(d_SDMPtrVect[i]->score);
          }
        }
      }
      // sort the SDM matrix by increasing scores
      std::sort(d_SDMPtrVect.begin(), d_SDMPtrVect.end(), this->compareSDMScore);
    }

    
    void SDM::fillFromDist(double threshold,const boost::dynamic_bitset<> &refHvyAtoms,
                           const boost::dynamic_bitset<> &prbHvyAtoms)
    {
      unsigned int n = 0;
      unsigned int pairs = 0;
      const RDGeom::POINT3D_VECT &refPos = d_refConf->getPositions();
      const RDGeom::POINT3D_VECT &prbPos = d_prbConf->getPositions();
      unsigned int mapIdx = 0;
      
      
      d_SDMPtrVect.clear();
      double sqThreshold = square(threshold);
      unsigned int refNAtoms = d_refConf->getNumAtoms();
      unsigned int prbNAtoms = d_prbConf->getNumAtoms();
      unsigned int largestNAtoms = std::max(prbNAtoms, refNAtoms);
      boost::dynamic_bitset<> refUsed(largestNAtoms);
      boost::dynamic_bitset<> prbUsed(largestNAtoms);
      // loop over ref atoms
      for (unsigned int i = 0; i < refNAtoms; ++i) {
        if (!refHvyAtoms[i]) continue;
        // loop over prb atoms
        for (unsigned int j = 0; j < prbNAtoms; ++j) {
          if (!prbHvyAtoms[j]) continue;
          double sqDist = (refPos[i] - prbPos[j]).lengthSq();
          // if the distance between these two atoms is lower
          // than threshold, then include this pair in the SDM matrix
          bool isConstraint = (d_o3aConstraintVect
            && (mapIdx < (*d_o3aConstraintVect).size())
            && (j == (*d_o3aConstraintVect)[mapIdx]->getPrbIdx())
            && (i == (*d_o3aConstraintVect)[mapIdx]->getRefIdx()));
          if ((sqDist < sqThreshold) || isConstraint) {
            SDMElement *sdmElement = new SDMElement();
            sdmElement->idx[0] = i;
            sdmElement->idx[1] = j;
            sdmElement->sqDist = sqDist;
            sdmElement->o3aConstraint = NULL;
            if (isConstraint) {
              sdmElement->o3aConstraint = (*d_o3aConstraintVect)[mapIdx];
              //in case there are duplicate constraints referring to the same
              //atom pair, apply only the one with the highest weight
              while ((mapIdx < (*d_o3aConstraintVect).size())
                && (j == (*d_o3aConstraintVect)[mapIdx]->getPrbIdx())
                && (i == (*d_o3aConstraintVect)[mapIdx]->getRefIdx())) {
                ++mapIdx;
              }
            }
            d_SDMPtrVect.push_back(boost::shared_ptr<SDMElement>(sdmElement));
            ++n;
          }
        }
      }
      // sort the SDM matrix by decreasing constraint weights and,
      // as a second criterion, increasing distances
      std::sort(d_SDMPtrVect.begin(), d_SDMPtrVect.end(), this->compareSDMDist);
      // increase the number of pairs which will be used
      // by rmsAlgorithm until an atom which has been
      // included in a previous pair is found
      for (unsigned int i = 0; (i < n) && (!(refUsed[d_SDMPtrVect[i]->idx[0]]))
        && (!(prbUsed[d_SDMPtrVect[i]->idx[1]])); ++i) {
        ++pairs;
        refUsed[d_SDMPtrVect[i]->idx[0]] = 1;
        prbUsed[d_SDMPtrVect[i]->idx[1]] = 1;
      }
      d_SDMPtrVect.resize(pairs);
    }
    

    double o3aMMFFWeightFunc(const unsigned int prbIdx, const unsigned int refIdx, void *data)
    {
      boost::uint8_t mmffSim = mmffSimMatrix
        [(static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFAtomType(refIdx) - 1]
        [(static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFAtomType(prbIdx) - 1];
      
      return static_cast<double>(O3_MAX_WEIGHT_COEFF
        - (static_cast<O3AFuncData *>(data))->weight) * ((1.0 + O3_CHARGE_COEFF
        * fabs((static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFPartialCharge(refIdx)
        + (static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFPartialCharge(prbIdx)))
        / (1.0 + fabs((static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFPartialCharge(refIdx)
        - (static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFPartialCharge(prbIdx))))
        + (double)((static_cast<O3AFuncData *>(data))->weight) / static_cast<double>(mmffSim);
    }


    double o3aCrippenWeightFunc(const unsigned int prbIdx, const unsigned int refIdx, void *data)
    {
      return (0.01 + fabs((*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->refProp)))[refIdx]
        + (*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->prbProp)))[prbIdx])
        / (1.0 + fabs((*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->refProp)))[refIdx]
        - (*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->prbProp)))[prbIdx])));
    }


    void SDM::prepareMatchWeightsVect(RDKit::MatchVectType &matchVect,
      RDNumeric::DoubleVector &weights, double (*weightFunc)
      (const unsigned int, const unsigned int, void *), void *data) {
      PRECONDITION(matchVect.size() == weights.size(), "matchVect/weights size mismatch");
      double min;
      for (unsigned int i = 0; i < matchVect.size(); ++i) {
        // first pair element is prb, second is ref
        matchVect[i].first = d_SDMPtrVect[i]->idx[1];
        matchVect[i].second = d_SDMPtrVect[i]->idx[0];
        weights[i] = (weightFunc ? weightFunc(d_SDMPtrVect[i]->idx[1],
          d_SDMPtrVect[i]->idx[0], data) : 1.0);
        if (d_SDMPtrVect[i]->o3aConstraint) {
          weights[i] += d_SDMPtrVect[i]->o3aConstraint->getWeight();
        }
        if ((!i) || (weights[i] < min)) {
          min = weights[i];
        }
      }
      for (unsigned int i = 0; i < weights.size(); ++i) {
        weights[i] /= min;
      }
    }
    

    double o3aMMFFScoringFunc(const unsigned int prbIdx,
      const unsigned int refIdx, void *data)
    {
      const RDGeom::POINT3D_VECT &refPos =
        (static_cast<O3AFuncData *>(data))->refConf->getPositions();
      const RDGeom::POINT3D_VECT &prbPos =
        (static_cast<O3AFuncData *>(data))->prbConf->getPositions();

      double sqDist = (refPos[refIdx] - prbPos[prbIdx]).lengthSq();
      double chargeDiff = fabs((static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFPartialCharge(refIdx)
        - (static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFPartialCharge(prbIdx));
      double chargeSum = fabs((static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->refProp))->getMMFFPartialCharge(refIdx)
        + (static_cast<MMFF::MMFFMolProperties *>
        ((static_cast<O3AFuncData *>(data))->prbProp))->getMMFFPartialCharge(prbIdx));

      return ((O3_SCORING_FUNCTION_ALPHA + (1.0 + O3_CHARGE_COEFF * chargeSum)
        / (1.0 + chargeDiff)) * exp(-O3_SCORING_FUNCTION_BETA * sqDist));
    }


    double o3aCrippenScoringFunc(const unsigned int prbIdx,
      const unsigned int refIdx, void *data)
    {
      const RDGeom::POINT3D_VECT &refPos =
        (static_cast<O3AFuncData *>(data))->refConf->getPositions();
      const RDGeom::POINT3D_VECT &prbPos =
        (static_cast<O3AFuncData *>(data))->prbConf->getPositions();

      double sqDist = (refPos[refIdx] - prbPos[prbIdx]).lengthSq();
      double logpDiff = fabs((*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->refProp)))[refIdx]
        - (*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->prbProp)))[prbIdx]);
      double logpSum = fabs((*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->refProp)))[refIdx]
        + (*(static_cast<std::vector<double> *>
        ((static_cast<O3AFuncData *>(data))->prbProp)))[prbIdx]);

      return ((O3_SCORING_FUNCTION_ALPHA + (1.0 + O3_CRIPPEN_COEFF * logpSum)
        / (1.0 + logpDiff)) * exp(-O3_SCORING_FUNCTION_BETA * sqDist));
    }


    double SDM::scoreAlignment(double (*scoringFunc)
      (const unsigned int, const unsigned int, void *), void *data)
    {
      double score = 0.0;
      
      for (unsigned int i = 0; i < d_SDMPtrVect.size(); ++i) {
        score += scoringFunc(d_SDMPtrVect[i]->idx[1], d_SDMPtrVect[i]->idx[0], data);
      }
      
      return score;
    }
    
    
    O3A::O3A(int (*costFunc)(const unsigned int, const unsigned int, double, void *),
      double (*weightFunc)(const unsigned int, const unsigned int, void *),
      double (*scoringFunc)(const unsigned int, const unsigned int, void *),
      void *data, ROMol &prbMol, const ROMol &refMol,
      const int prbCid, const int refCid,
      boost::dynamic_bitset<> *prbHvyAtoms,
      boost::dynamic_bitset<> *refHvyAtoms,
      const bool reflect, const unsigned int maxIters, unsigned int options,
      O3AConstraintVect *o3aConstraintVect, ROMol *extWorkPrbMol,
      LAP *extLAP, MolHistogram *extPrbHist, MolHistogram *extRefHist) :
      d_prbMol(&prbMol),
      d_refMol(&refMol),
      d_prbCid(prbCid),
      d_refCid(refCid),
      d_reflect(reflect),
      d_maxIters(maxIters)
    {
      ROMol *workPrbMol = (extWorkPrbMol ? extWorkPrbMol : new ROMol(prbMol));
      Conformer &prbConf = workPrbMol->getConformer();
      const Conformer &refConf = refMol.getConformer();
      unsigned int accuracy = options & O3_ACCURACY_MASK;
      bool local = options & O3_LOCAL_ONLY;
      unsigned int refNAtoms = refMol.getNumAtoms();
      unsigned int prbNAtoms = prbMol.getNumAtoms();
      unsigned int refNHeavyAtoms = refMol.getNumHeavyAtoms();
      unsigned int prbNHeavyAtoms = prbMol.getNumHeavyAtoms();
      unsigned int largestNHeavyAtoms = std::max(refNHeavyAtoms, prbNHeavyAtoms);
      unsigned int i;
      if (!refHvyAtoms) {
        refHvyAtoms = new boost::dynamic_bitset<>(refNAtoms);
        for (i = 0; i < refNAtoms; ++i) {
          if (refMol[i]->getAtomicNum() != 1) {
            refHvyAtoms->set(i);
          }
        }
      }
      if (!prbHvyAtoms) {
        for (i = 0; i < prbNAtoms; ++i) {
          if (prbMol[i]->getAtomicNum() != 1) {
            prbHvyAtoms->set(i);
          }
        }
      }
      std::vector<unsigned int> pairs(4, 0);
      std::vector<double> score(3, 0.0);
      std::vector<double> pairsRMSD(2, 0.0);
      std::vector<SDM *> bestSDM((unsigned int)3, NULL);
      SDM startSDM(&prbConf, &refConf, o3aConstraintVect);
      MolHistogram *refHist = NULL;
      MolHistogram *prbHist = NULL;
      LAP *lap = NULL;
      if (local) {
        startSDM.fillFromDist(O3_SDM_THRESHOLD_START, *refHvyAtoms, *prbHvyAtoms);
      }
      else {
        refHist = (extRefHist ? extRefHist
          : new MolHistogram(refMol, MolOps::get3DDistanceMat(refMol, refCid)));
        prbHist = (extPrbHist ? extPrbHist
          : new MolHistogram(prbMol, MolOps::get3DDistanceMat(prbMol, prbCid)));
        lap = (extLAP ? extLAP : new LAP(largestNHeavyAtoms));
        lap->computeCostMatrix(prbMol, *prbHist,
          refMol, *refHist, o3aConstraintVect, costFunc, data);
        lap->computeMinCostPath(largestNHeavyAtoms);
        startSDM.fillFromLAP(*lap);
      }
      for (pairs[0] = 0, score[0] = 0.0, i = 3; i < startSDM.size(); ++i) {
        RDKit::MatchVectType startMatchVect(i);
        RDNumeric::DoubleVector startWeights(i);
        startSDM.prepareMatchWeightsVect(startMatchVect,
          startWeights, weightFunc, data);
        alignMol(*workPrbMol, refMol, prbCid, refCid,
          &startMatchVect, &startWeights, reflect, maxIters);
        unsigned int sdmThresholdIt;
        for (sdmThresholdIt = ((accuracy < 1) ? 0 : O3_MAX_SDM_THRESHOLD_ITER - 1),
          pairs[1] = 0, score[1] = 0.0; sdmThresholdIt < O3_MAX_SDM_THRESHOLD_ITER;
          ++sdmThresholdIt) {
          pairs[2] = 0;
          bool flag = true;
          unsigned int iter = 0;
          double sdmThresholdDist = O3_SDM_THRESHOLD_START
            + (double)sdmThresholdIt * O3_SDM_THRESHOLD_STEP;
          while (flag && (iter < O3_MAX_SDM_ITERATIONS)) {
            SDM progressSDM(&prbConf, &refConf, o3aConstraintVect);
            progressSDM.fillFromDist(sdmThresholdDist, *refHvyAtoms, *prbHvyAtoms);
            pairs[3] = progressSDM.size();
            if (pairs[3] < 3) {
              break;
            }
            RDKit::MatchVectType progressMatchVect(pairs[3]);
            RDNumeric::DoubleVector progressWeights(pairs[3]);
            progressSDM.prepareMatchWeightsVect(progressMatchVect,
              progressWeights, weightFunc, data);
            RDGeom::Transform3D trans;
            pairsRMSD[1] = getAlignmentTransform(*workPrbMol, refMol,
              trans, prbCid, refCid, &progressMatchVect, &progressWeights,
              reflect, maxIters);
            flag = ((pairs[3] > pairs[2]) || ((pairs[3] == pairs[2])
              && ((!iter) || ((pairsRMSD[0] - pairsRMSD[1]) > O3_RMSD_THRESHOLD))));
            if (flag) {
              pairs[2] = pairs[3];
              pairsRMSD[0] = pairsRMSD[1];
              if (bestSDM[2]) {
                delete bestSDM[2];
              }
              bestSDM[2] = new SDM();
              *(bestSDM[2]) = progressSDM;
              MolTransforms::transformConformer(prbConf, trans);
            }
            ++iter;
          }
          score[2] = ((pairs[2] >= 3)
            ? bestSDM[2]->scoreAlignment(scoringFunc, data) : 0.0);
          if ((score[2] - score[1]) > O3_SCORE_THRESHOLD) {
            pairs[1] = pairs[2];
            score[1] = score[2];
            if (bestSDM[1]) {
              delete bestSDM[1];
            }
            bestSDM[1] = new SDM();
            *(bestSDM[1]) = *(bestSDM[2]);
          }
        }
        if ((score[1] - score[0]) > O3_SCORE_THRESHOLD) {
          pairs[0] = pairs[1];
          score[0] = score[1];
          if (bestSDM[0]) {
            delete bestSDM[0];
          }
          bestSDM[0] = new SDM();
          *(bestSDM[0]) = *(bestSDM[1]);
        }
      }
      if ((pairs[0] < 3) && (startSDM.size() >= 3)) {
        pairs[0] = startSDM.size();
        if (bestSDM[0]) {
          delete bestSDM[0];
        }
        bestSDM[0] = new SDM();
        *(bestSDM[0]) = startSDM;
        score[0] = bestSDM[0]->scoreAlignment(scoringFunc, data);
      }
      RDKit::MatchVectType *o3aMatchVect = new RDKit::MatchVectType(pairs[0]);
      RDNumeric::DoubleVector *o3aWeights = new RDNumeric::DoubleVector(pairs[0]);
      d_o3aMatchVect = o3aMatchVect;
      d_o3aWeights = o3aWeights;
      if (pairs[0] >= 3) {
        bestSDM[0]->prepareMatchWeightsVect(*o3aMatchVect, *o3aWeights, weightFunc, data);
        d_o3aScore = score[0];
      }
      else {
        d_o3aScore = 0.0;
      }
      for (i = 0; i < 3; ++i) {
        if (bestSDM[i]) {
          delete bestSDM[i];
        }
      }
      if ((!extLAP) && lap) {
        delete lap;
      }
      if ((!extRefHist) && refHist) {
        delete refHist;
      }
      if ((!extPrbHist) && prbHist) {
        delete prbHist;
      }
      if ((!extWorkPrbMol) && workPrbMol) {
        delete workPrbMol;
      }
    }


    O3A::O3A(ROMol &prbMol, const ROMol &refMol, void *prbProp, void *refProp,
      AtomTypeScheme atomTypes, const int prbCid, const int refCid,
      const bool reflect, const unsigned int maxIters, unsigned int options,
      const MatchVectType *constraintMap, const RDNumeric::DoubleVector *constraintWeights,
      LAP *extLAP, MolHistogram *extPrbHist, MolHistogram *extRefHist) :
      d_prbMol(&prbMol),
      d_refMol(&refMol),
      d_prbCid(prbCid),
      d_refCid(refCid),
      d_reflect(reflect),
      d_maxIters(maxIters)
    {
      int (*costFunc)(const unsigned int, const unsigned int, double, void *) = o3aMMFFCostFunc;
      double (*weightFunc)(const unsigned int, const unsigned int, void *) = o3aMMFFWeightFunc;
      double (*scoringFunc)(const unsigned int, const unsigned int, void *) = o3aMMFFScoringFunc;
      if (atomTypes == CRIPPEN) {
        costFunc = o3aCrippenCostFunc;
        weightFunc = o3aCrippenWeightFunc;
        scoringFunc = o3aCrippenScoringFunc;
      }
      O3AFuncData data;
      ROMol extWorkPrbMol(prbMol);
      data.prbConf = &(extWorkPrbMol.getConformer());
      data.refConf = &(refMol.getConformer());
      data.prbProp = prbProp;
      data.refProp = refProp;
      unsigned int refNAtoms = refMol.getNumAtoms();
      unsigned int prbNAtoms = prbMol.getNumAtoms();
      boost::dynamic_bitset<> refHvyAtoms(refNAtoms);
      boost::dynamic_bitset<> prbHvyAtoms(prbNAtoms);
      unsigned int i;
      for (i = 0; i < refNAtoms; ++i) {
        if (refMol[i]->getAtomicNum() != 1) {
          refHvyAtoms.set(i);
        }
      }
      for (i = 0; i < prbNAtoms; ++i) {
        if (prbMol[i]->getAtomicNum() != 1) {
          prbHvyAtoms.set(i);
        }
      }
      unsigned int refNHeavyAtoms = refMol.getNumHeavyAtoms();
      unsigned int prbNHeavyAtoms = prbMol.getNumHeavyAtoms();
      unsigned int largestNHeavyAtoms = std::max(refNHeavyAtoms, prbNHeavyAtoms);
      unsigned int accuracy = options & O3_ACCURACY_MASK;
      bool local = options & O3_LOCAL_ONLY;
      O3AConstraintVect o3aConstraintVect;
      if (constraintMap) {
        if (constraintWeights) {
          if ((*constraintMap).size() != (*constraintWeights).size()) {
            throw MolAlignException("The number of weights should match the number of constraints");
          }
        }
        for (i = 0; i < (*constraintMap).size(); ++i) {
          if (((*constraintMap)[i].first < 0) || ((*constraintMap)[i].first >= prbNAtoms)
            || ((*constraintMap)[i].second < 0) || ((*constraintMap)[i].second >= refNAtoms)) {
            throw MolAlignException("Constrained atom idx out of range");
          }
          if ((!prbHvyAtoms[(*constraintMap)[i].first]) || (!refHvyAtoms[(*constraintMap)[i].second])) {
            throw MolAlignException("Constrained atoms must be heavy atoms");
          }
          o3aConstraintVect.append((*constraintMap)[i].first, (*constraintMap)[i].second,
            constraintWeights ? (*constraintWeights)[i] : O3_DEFAULT_CONSTRAINT_WEIGHT);
        }
      }
      int c;
      int l;
      std::vector<double> score(2, 0.0);
      O3A *bestO3A = NULL;
      MolHistogram *refHist = NULL;
      MolHistogram *prbHist = NULL;
      LAP *lap = NULL;
      if (!local) {
        refHist = (extRefHist ? extRefHist
          : new MolHistogram(refMol, MolOps::get3DDistanceMat(refMol, refCid)));
        prbHist = (extPrbHist ? extPrbHist
          : new MolHistogram(prbMol, MolOps::get3DDistanceMat(prbMol, prbCid)));
        lap = (extLAP ? extLAP : new LAP(largestNHeavyAtoms));
      }
      for (l = 0, score[0] = 0.0;
        l <= (((accuracy < 2) && (atomTypes == MMFF94)) ? 1 : 0); ++l) {
        data.useMMFFSim = (bool)l;
        for (c = 0; c <= O3_MAX_WEIGHT_COEFF;
          c += ((accuracy < 3) ? 1 : (O3_MAX_WEIGHT_COEFF + 1))) {
          data.weight = (l ? c : 0);
          data.coeff = c;
          O3A *o3a = new O3A(costFunc, weightFunc, scoringFunc,
            &data, prbMol, refMol, prbCid, refCid, &prbHvyAtoms, &refHvyAtoms,
            reflect, maxIters, options, &o3aConstraintVect, &extWorkPrbMol,
            lap, prbHist, refHist);
          score[1] = o3a->score();
          if (((score[1] - score[0]) > O3_SCORE_THRESHOLD)
              || isDoubleZero(score[0])) {
            if (bestO3A) {
              delete bestO3A;
            }
            bestO3A = o3a;
            score[0] = score[1];
          }
          else {
            delete o3a;
          }
        }
      }
      unsigned int pairs = bestO3A->matches()->size();
      RDKit::MatchVectType *bestO3AMatchVect = new RDKit::MatchVectType(pairs);
      RDNumeric::DoubleVector *bestO3AWeights = new RDNumeric::DoubleVector(pairs);
      d_o3aMatchVect = bestO3AMatchVect;
      d_o3aWeights = bestO3AWeights;
      if (pairs >= 3) {
        for (i = 0; i < pairs; ++i) {
          (*bestO3AMatchVect)[i].first = (*(bestO3A->matches()))[i].first;
          (*bestO3AMatchVect)[i].second = (*(bestO3A->matches()))[i].second;
          (*bestO3AWeights)[i] = (*(bestO3A->weights()))[i];
        }
        d_o3aScore = bestO3A->score();
      }
      else {
        d_o3aScore = 0.0;
      }
      delete bestO3A;
      if (!extLAP) {
        delete lap;
      }
      if (!extRefHist) {
        delete refHist;
      }
      if (!extPrbHist) {
        delete prbHist;
      }
    }


    double _rmsdMatchVect(ROMol *d_prbMol, const ROMol *d_refMol, 
      const RDGeom::POINT3D_VECT &prbPos, const RDGeom::POINT3D_VECT &refPos,
      const RDKit::MatchVectType *matchVect)
    {
      double rmsd = 0.0;
      if (matchVect) {
        for (unsigned int i = 0; i < (*matchVect).size(); ++i) {
          //first pair element is prb, second is ref
          rmsd += (prbPos[(*matchVect)[i].first]
            - refPos[(*matchVect)[i].second]).lengthSq();
        }
        rmsd = sqrt(rmsd / (*matchVect).size());
      }
      else {
        RDGeom::Point3D refCtd;
        RDGeom::Point3D prbCtd;
        unsigned int i;
        unsigned int j;
        unsigned int nHeavy;
        for (i = 0, nHeavy = 0; i < refPos.size(); ++i) {
          if (d_refMol->getAtomWithIdx(i)->getAtomicNum() == 1) {
            continue;
          }
          refCtd[0] += refPos[i][0];
          refCtd[1] += refPos[i][1];
          refCtd[2] += refPos[i][2];
          ++nHeavy;
        }
        refCtd[0] /= (double)nHeavy;
        refCtd[1] /= (double)nHeavy;
        refCtd[2] /= (double)nHeavy;
        for (i = 0, nHeavy = 0; i < prbPos.size(); ++i) {
          if (d_prbMol->getAtomWithIdx(i)->getAtomicNum() == 1) {
            continue;
          }
          prbCtd[0] += prbPos[i][0];
          prbCtd[1] += prbPos[i][1];
          prbCtd[2] += prbPos[i][2];
          ++nHeavy;
        }
        prbCtd[0] /= (double)nHeavy;
        prbCtd[1] /= (double)nHeavy;
        prbCtd[2] /= (double)nHeavy;
        rmsd = (refCtd - prbCtd).length();
      }
      
      return rmsd;
    };


    double O3A::align() {
      double rmsd = 0.0;
      const RDKit::MatchVectType *matchVectPtr = NULL;
      
      if (d_o3aMatchVect && (d_o3aMatchVect->size() >= 3)) {
        alignMol(*d_prbMol, *d_refMol, d_prbCid, d_refCid,
          d_o3aMatchVect, d_o3aWeights, d_reflect, d_maxIters);
        matchVectPtr = d_o3aMatchVect;
      }
      rmsd = _rmsdMatchVect(d_prbMol, d_refMol, 
        d_prbMol->getConformer(d_prbCid).getPositions(),
        d_refMol->getConformer(d_refCid).getPositions(), matchVectPtr);

      return rmsd;
    }


    double O3A::trans(RDGeom::Transform3D &trans) {
      double rmsd = 0.0;
      const RDKit::MatchVectType *matchVectPtr = NULL;
      Conformer transConf(d_prbMol->getConformer(d_prbCid));
      if (d_o3aMatchVect && (d_o3aMatchVect->size() >= 3)) {
        getAlignmentTransform(*d_prbMol, *d_refMol,
          trans, d_prbCid, d_refCid, d_o3aMatchVect, d_o3aWeights,
          d_reflect, d_maxIters);
        MolTransforms::transformConformer(transConf, trans);
        matchVectPtr = d_o3aMatchVect;
      }
      rmsd = _rmsdMatchVect(d_prbMol, d_refMol, transConf.getPositions(),
        d_refMol->getConformer(d_refCid).getPositions(), d_o3aMatchVect);

      return rmsd;
    };


    void randomTransform(ROMol &mol, const int cid, const int seed)
    {
      if (seed > 0) {
        rng_type &generator = getRandomGenerator();
        generator.seed(seed);
      }
      Conformer &conf = mol.getConformer(cid);
      RDGeom::POINT3D_VECT &pos = conf.getPositions();
      double maxCoord[3];
      double minCoord[3];
      RDGeom::Point3D rndTrans;
      double rndRot[3];
      for (unsigned int j = 0; j < 3; ++j) {
        for (unsigned int i = 0; i < pos.size(); ++i) {
          if ((!i) || (pos[i][j] > maxCoord[j])) {
            maxCoord[j] = pos[i][j];
          }
          if ((!i) || (pos[i][j] < minCoord[j])) {
            minCoord[j] = pos[i][j];
          }
        }
        rndTrans[j] = (maxCoord[j] - minCoord[j])
          * O3_RANDOM_TRANS_COEFF * (getRandomVal() - 0.5);
        rndRot[j] = getRandomVal() * 2.0 * M_PI;
      }
      RDGeom::Point3D ctd = MolTransforms::computeCentroid(conf, false);
      RDGeom::Transform3D transToOrigin;
      RDGeom::Transform3D xRot;
      RDGeom::Transform3D yRot;
      RDGeom::Transform3D zRotAndTrans;
      transToOrigin.SetTranslation(-ctd);
      xRot.SetRotation(rndRot[0], RDGeom::X_Axis);
      yRot.SetRotation(rndRot[1], RDGeom::Y_Axis);
      zRotAndTrans.SetRotation(rndRot[2], RDGeom::Z_Axis);
      zRotAndTrans.SetTranslation(ctd + rndTrans);
      MolTransforms::transformConformer(conf,
        zRotAndTrans * yRot * xRot * transToOrigin);
    }
  }
}