//
// Copyright (C) 2013-2018 Paolo Tosco, Greg Landrum
//
// @@ 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>
#include <RDGeneral/RDThreads.h>

#ifdef RDK_THREADSAFE_SSS
#include <thread>
#include <future>
#endif

#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,
                           bool cleanupDmat)
    : 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;
  }
  if (cleanupDmat) delete[] dmat;
}

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;
    }
    auto *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 = nullptr;
    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) {
        auto *sdmElement = new SDMElement();
        sdmElement->idx[0] = i;
        sdmElement->idx[1] = j;
        sdmElement->sqDist = sqDist;
        sdmElement->o3aConstraint = nullptr;
        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 = MAX_DOUBLE;
  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 (auto &ptr : d_SDMPtrVect) {
    score += scoringFunc(ptr->idx[1], ptr->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(prbCid);
  const Conformer &refConf = refMol.getConformer(refCid);
  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, nullptr);
  SDM startSDM(&prbConf, &refConf, o3aConstraintVect);
  MolHistogram *refHist = nullptr;
  MolHistogram *prbHist = nullptr;
  LAP *lap = nullptr;
  if (local) {
    startSDM.fillFromDist(O3_SDM_THRESHOLD_START, *refHvyAtoms, *prbHvyAtoms);
  } else {
    refHist = (extRefHist ? extRefHist
                          : new MolHistogram(
                                refMol, MolOps::get3DDistanceMat(
                                            refMol, refCid, false, false, ""),
                                true));
    prbHist = (extPrbHist ? extPrbHist
                          : new MolHistogram(
                                prbMol, MolOps::get3DDistanceMat(
                                            prbMol, prbCid, false, false, ""),
                                true));
    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]);
  auto *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(prbCid));
  data.refConf = &(refMol.getConformer(refCid));
  data.prbProp = prbProp;
  data.refProp = refProp;
  int refNAtoms = rdcast<int>(refMol.getNumAtoms());
  int prbNAtoms = rdcast<int>(prbMol.getNumAtoms());
  boost::dynamic_bitset<> refHvyAtoms(refNAtoms);
  boost::dynamic_bitset<> prbHvyAtoms(prbNAtoms);

  for (int i = 0; i < refNAtoms; ++i) {
    if (refMol[i]->getAtomicNum() != 1) {
      refHvyAtoms.set(i);
    }
  }
  for (int 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 (unsigned int 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 = nullptr;
  MolHistogram *refHist = nullptr;
  MolHistogram *prbHist = nullptr;
  LAP *lap = nullptr;
  if (!local) {
    refHist = (extRefHist ? extRefHist
                          : new MolHistogram(
                                refMol, MolOps::get3DDistanceMat(
                                            refMol, refCid, false, false, ""),
                                true));
    prbHist = (extPrbHist ? extPrbHist
                          : new MolHistogram(
                                prbMol, MolOps::get3DDistanceMat(
                                            prbMol, prbCid, false, false, ""),
                                true));
    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;
      auto *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 = rdcast<unsigned int>(bestO3A->matches()->size());
  RDKit::MatchVectType *bestO3AMatchVect = new RDKit::MatchVectType(pairs);
  auto *bestO3AWeights = new RDNumeric::DoubleVector(pairs);
  d_o3aMatchVect = bestO3AMatchVect;
  d_o3aWeights = bestO3AWeights;
  if (pairs >= 3) {
    for (unsigned int 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 (const auto &i : (*matchVect)) {
      // first pair element is prb, second is ref
      rmsd += (prbPos[i.first] - refPos[i.second]).lengthSq();
    }
    rmsd = sqrt(rmsd / (*matchVect).size());
  } else {
    RDGeom::Point3D refCtd;
    RDGeom::Point3D prbCtd;
    unsigned int i;
    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 = nullptr;

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

#ifdef RDK_THREADSAFE_SSS
namespace detail {
typedef struct {
  O3A::AtomTypeScheme atomTypes;
  int refCid;
  bool reflect;
  unsigned int maxIters;
  unsigned int options;
  MatchVectType const *constraintMap;
  RDNumeric::DoubleVector const *constraintWeights;
} O3AHelperArgs_;
void O3AHelper_(ROMol *prbMol, const ROMol *refMol, void *prbProp,
                void *refProp, std::vector<boost::shared_ptr<O3A>> *res,
                unsigned int threadIdx, int numThreads,
                const O3AHelperArgs_ *args) {
  unsigned int i = 0;
  for (ROMol::ConstConformerIterator cit = prbMol->beginConformers();
       cit != prbMol->endConformers(); ++cit, ++i) {
    if (i % numThreads != threadIdx) continue;

    auto *lres =
        new O3A(*prbMol, *refMol, prbProp, refProp, args->atomTypes,
                (*cit)->getId(), args->refCid, args->reflect, args->maxIters,
                args->options, args->constraintMap, args->constraintWeights);
    (*res)[i].reset(lres);
  }
}
}  // end of detail namespace
#endif

void getO3AForProbeConfs(ROMol &prbMol, const ROMol &refMol, void *prbProp,
                         void *refProp,
                         std::vector<boost::shared_ptr<O3A>> &res,
                         int numThreads, O3A::AtomTypeScheme atomTypes,
                         const int refCid, const bool reflect,
                         const unsigned int maxIters, unsigned int options,
                         const MatchVectType *constraintMap,
                         const RDNumeric::DoubleVector *constraintWeights) {
  numThreads = getNumThreadsToUse(numThreads);
  res.resize(prbMol.getNumConformers());
  if (numThreads == 1) {
    unsigned int i = 0;
    for (ROMol::ConstConformerIterator cit = prbMol.beginConformers();
         cit != prbMol.endConformers(); ++cit, ++i) {
      auto *lres = new O3A(prbMol, refMol, prbProp, refProp, atomTypes,
                           (*cit)->getId(), refCid, reflect, maxIters, options,
                           constraintMap, constraintWeights);
      res[i].reset(lres);
    }
  }
#ifdef RDK_THREADSAFE_SSS
  else {
    std::vector<std::future<void>> tg;
    detail::O3AHelperArgs_ args = {atomTypes,        refCid,  reflect,
                                   maxIters,         options, constraintMap,
                                   constraintWeights};
    for (int ti = 0; ti < numThreads; ++ti) {
      tg.emplace_back(std::async(std::launch::async, detail::O3AHelper_,
                                 &prbMol, &refMol, prbProp, refProp, &res, ti,
                                 numThreads, &args));
    }
    for (auto &fut : tg) {
      fut.get();
    }
  }
#endif
}

}  // end of namespace MolAlign
}  // end of namespace RDKit