sketcherMinimizerFragment.cpp

/*
 Contributors: Nicola Zonta
 Copyright Schrodinger, LLC. All rights reserved
 */

#include "sketcherMinimizerFragment.h"
#include "sketcherMinimizerAtom.h"
#include "sketcherMinimizerBond.h"
#include "sketcherMinimizerMaths.h"
#include "sketcherMinimizerRing.h"

static const float INVERTED_MACROCYCLE_BOND_PENALTY = 100.f;
static const float SCALE_FRAGMENT_PENALTY = 500.f;
static const float SCALE_ATOMS_PENALTY = 50.f;
static const float ROTATE_FRAGMENT_PENALTY = 400.f;
static const float BREAK_CHAIN_PENALTY = 10.f;
static const float CHANGE_PARENT_BOND_PENALTY = 200.f;
static const float FLIP_RING_PENALTY = 200.f;
static const float FLIP_CONSTRAINED_FRAGMENT_PENALTY = 1000.f;

static const int FLIP_FRAGMENT_TIER = 0;
static const int INVERT_BOND_TIER = 1;
static const int ROTATE_FRAGMENT_TIER = 3;
static const int ELONGATE_PARENT_BOND_TIER = 2;
static const int SCALE_ATOMS_TIER = 4;
static const int FLIP_RING_TIER = 1;
static const int SCALE_FRAGMENT_TIER = 5;

CoordgenFragmentDOF::CoordgenFragmentDOF(sketcherMinimizerFragment* fragment)
    : m_currentState(0), m_optimalState(0), m_fragment(fragment)
{
}

CoordgenFragmentDOF::~CoordgenFragmentDOF() = default;

short unsigned int CoordgenFragmentDOF::getCurrentState()
{
    return m_currentState;
}

void CoordgenFragmentDOF::setState(short unsigned int state)
{
    m_currentState = state;
}

void CoordgenFragmentDOF::storeCurrentValueAsOptimal()
{
    m_optimalState = m_currentState;
}

void CoordgenFragmentDOF::setToOptimalValue()
{
    m_currentState = m_optimalState;
}

void CoordgenFragmentDOF::changeState()
{
    m_currentState++;
    m_currentState = m_currentState % numberOfStates();
}

void CoordgenFragmentDOF::addAtom(sketcherMinimizerAtom* atom)
{
    m_atoms.push_back(atom);
    atom->fragment->addDofToAtom(atom, this);
}

float CoordgenFragmentDOF::getCurrentPenalty() const
{
    return 0.f;
}

sketcherMinimizerFragment* CoordgenFragmentDOF::getFragment() const
{
    return m_fragment;
}

CoordgenFlipFragmentDOF::CoordgenFlipFragmentDOF(
    sketcherMinimizerFragment* fragment)
    : CoordgenFragmentDOF(fragment)
{
}

float CoordgenFlipFragmentDOF::getCurrentPenalty() const
{
    float penalty = 0.f;
    if (m_currentState != 0 && m_fragment->constrainedFlip) {
        penalty += FLIP_CONSTRAINED_FRAGMENT_PENALTY;
    }
    if (m_fragment->isChain && m_fragment->getParent() &&
        m_fragment->getParent()->isChain) {
        penalty += BREAK_CHAIN_PENALTY;
    }
    return penalty;
}

int CoordgenFlipFragmentDOF::numberOfStates() const
{
    if (m_fragment->getParent() == nullptr) {
        return 1;
    }
    return 2;
}

int CoordgenFlipFragmentDOF::tier() const
{
    return FLIP_FRAGMENT_TIER;
}

void CoordgenFlipFragmentDOF::apply() const
{
    if (m_currentState != 0) {
        for (auto& atom : m_fragment->_coordinates) {
            atom.first->coordinates.setY(-atom.first->coordinates.y());
        }
    }
}

CoordgenScaleFragmentDOF::CoordgenScaleFragmentDOF(
    sketcherMinimizerFragment* fragment)
    : CoordgenFragmentDOF(fragment)
{
}

int CoordgenScaleFragmentDOF::numberOfStates() const
{
    if (m_fragment->getRings().empty()) {
        return 1;
    }
    return 5;
}

int CoordgenScaleFragmentDOF::tier() const
{
    return SCALE_FRAGMENT_TIER;
}

void CoordgenScaleFragmentDOF::apply() const
{
    if (m_currentState != 0) {
        auto scale = static_cast<float>(pow(1.4, (m_currentState + 1) / 2));
        if (m_currentState % 2 == 0) {
            scale = 1 / scale;
        }
        for (auto& atom : m_fragment->_coordinates) {
            atom.first->setCoordinates(atom.first->getCoordinates() * scale);
        }
    }
}

float CoordgenScaleFragmentDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return SCALE_FRAGMENT_PENALTY * ((m_currentState + 1) / 2);
    }
    return 0.f;
}

CoordgenScaleAtomsDOF::CoordgenScaleAtomsDOF(sketcherMinimizerAtom* pivotAtom)
    : CoordgenFragmentDOF(pivotAtom->getFragment()), m_pivotAtom(pivotAtom)
{
}

int CoordgenScaleAtomsDOF::numberOfStates() const
{
    return 2;
}

int CoordgenScaleAtomsDOF::tier() const
{
    return SCALE_ATOMS_TIER;
}

void CoordgenScaleAtomsDOF::apply() const
{
    if (m_currentState != 0) {
        float scale = 0.4f;
        for (auto atom : m_atoms) {
            auto distance =
                atom->getCoordinates() - m_pivotAtom->getCoordinates();
            atom->setCoordinates(distance * scale +
                                 m_pivotAtom->getCoordinates());
        }
    }
}

float CoordgenScaleAtomsDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return SCALE_ATOMS_PENALTY * m_atoms.size();
    }
    return 0.f;
}

CoordgenChangeParentBondLengthFragmentDOF::
    CoordgenChangeParentBondLengthFragmentDOF(
        sketcherMinimizerFragment* fragment)
    : CoordgenFragmentDOF(fragment)
{
}

int CoordgenChangeParentBondLengthFragmentDOF::numberOfStates() const
{
    return 7;
}

int CoordgenChangeParentBondLengthFragmentDOF::tier() const
{
    return ELONGATE_PARENT_BOND_TIER;
}

void CoordgenChangeParentBondLengthFragmentDOF::apply() const
{
    if (m_currentState != 0) {
        auto scale = static_cast<float>(pow(1.6, (m_currentState + 1) / 2));
        if (m_currentState % 2 == 0) {
            scale = 1 / scale;
        }
        float MoveBy = BONDLENGTH * (scale - 1);
        for (auto& atom : m_fragment->_coordinates) {
            atom.first->coordinates.setX(atom.first->coordinates.x() + MoveBy);
        }
    }
}

float CoordgenChangeParentBondLengthFragmentDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return CHANGE_PARENT_BOND_PENALTY * ((m_currentState + 1) / 2);
    }
    return 0.f;
}

CoordgenRotateFragmentDOF::CoordgenRotateFragmentDOF(
    sketcherMinimizerFragment* fragment)
    : CoordgenFragmentDOF(fragment)
{
}

int CoordgenRotateFragmentDOF::numberOfStates() const
{
    if (m_fragment->getParent() == nullptr) {
        return 1;
    }
    return 5;
}

int CoordgenRotateFragmentDOF::tier() const
{
    return ROTATE_FRAGMENT_TIER;
}

void CoordgenRotateFragmentDOF::apply() const
{
    if (m_currentState != 0) {
        auto angle =
            static_cast<float>(M_PI / 180 * 15 * ((m_currentState + 1) / 2));
        if (m_currentState % 2 == 0) {
            angle = -angle;
        }
        float sine = sin(angle);
        float cosine = cos(angle);
        for (auto& atom : m_fragment->_coordinates) {
            sketcherMinimizerPointF origin(-BONDLENGTH, 0);
            sketcherMinimizerPointF coords =
                atom.first->getCoordinates() - origin;
            coords.rotate(sine, cosine);
            atom.first->setCoordinates(coords + origin);
        }
    }
}

float CoordgenRotateFragmentDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return ROTATE_FRAGMENT_PENALTY * ((m_currentState + 1) / 2);
    }
    return 0.f;
}

CoordgenInvertBondDOF::CoordgenInvertBondDOF(sketcherMinimizerAtom* pivotAtom,
                                             sketcherMinimizerAtom* boundAtom)
    : CoordgenFragmentDOF(pivotAtom->getFragment()), m_pivotAtom(pivotAtom),
      m_boundAtom(boundAtom)
{
    assert(pivotAtom->bondTo(boundAtom) != nullptr);
    addAtom(boundAtom);
}

int CoordgenInvertBondDOF::numberOfStates() const
{
    return 2;
}

int CoordgenInvertBondDOF::tier() const
{
    return INVERT_BOND_TIER;
}

void CoordgenInvertBondDOF::apply() const
{
    if (m_currentState != 0) {
        sketcherMinimizerPointF pivot = m_pivotAtom->getCoordinates();
        sketcherMinimizerPointF bondDir = m_boundAtom->getCoordinates() - pivot;
        sketcherMinimizerPointF normal(bondDir.y(), -bondDir.x());
        sketcherMinimizerPointF point1 = pivot + normal;
        sketcherMinimizerPointF point2 = pivot - normal;

        for (auto& atom : m_atoms) {
            atom->setCoordinates(sketcherMinimizerMaths::mirrorPoint(
                atom->getCoordinates(), point1, point2));
        }
    }
}

float CoordgenInvertBondDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return INVERTED_MACROCYCLE_BOND_PENALTY;
    }
    return 0.f;
}

CoordgenFlipRingDOF::CoordgenFlipRingDOF(
    sketcherMinimizerRing* ring,
    const std::vector<sketcherMinimizerAtom*>& fusionAtoms)
    : CoordgenFragmentDOF((*fusionAtoms.begin())->getFragment()),
      m_pivotAtom1(*(fusionAtoms.begin())),
      m_pivotAtom2(*(fusionAtoms.rbegin())),
      m_penalty(std::abs(int(ring->size() - 2 * fusionAtoms.size() + 2)))
{
    for (auto atom : ring->getAtoms()) {
        addAtom(atom);
    }
}

int CoordgenFlipRingDOF::numberOfStates() const
{
    return 2;
}

int CoordgenFlipRingDOF::tier() const
{
    return FLIP_RING_TIER;
}

void CoordgenFlipRingDOF::apply() const
{
    if (m_currentState != 0) {
        for (auto& atom : m_atoms) {
            atom->setCoordinates(sketcherMinimizerMaths::mirrorPoint(
                atom->getCoordinates(), m_pivotAtom1->getCoordinates(),
                m_pivotAtom2->getCoordinates()));
        }
    }
}

float CoordgenFlipRingDOF::getCurrentPenalty() const
{
    if (m_currentState != 0) {
        return FLIP_RING_PENALTY * m_penalty;
    }
    return 0.f;
}

sketcherMinimizerFragment::sketcherMinimizerFragment()
    : fixed(false), isTemplated(false), constrained(false), isChain(false),
      _bondToParent(nullptr), longestChainFromHere(0.f),
      numberOfChildrenAtoms(0), numberOfChildrenAtomsRank(0.f),
      m_parent(nullptr)
{
    m_dofs.push_back(new CoordgenFlipFragmentDOF(this));
    //    m_dofs.push_back(new CoordgenScaleFragmentDOF(this));
    m_dofs.push_back(new CoordgenChangeParentBondLengthFragmentDOF(this));
    m_dofs.push_back(new CoordgenRotateFragmentDOF(this));
}

sketcherMinimizerFragment::~sketcherMinimizerFragment()
{
    for (auto dof : m_dofs) {
        delete dof;
    }
}

void sketcherMinimizerFragment::addDof(CoordgenFragmentDOF* dof)
{
    m_dofs.push_back(dof);
}

const std::vector<CoordgenFragmentDOF*>& sketcherMinimizerFragment::getDofs()
{
    return m_dofs;
}

unsigned int sketcherMinimizerFragment::totalWeight() const
{
    int n = 0;
    for (auto m_atom : m_atoms) {
        n += m_atom->atomicNumber + m_atom->_implicitHs;
    }
    return n;
}
unsigned int sketcherMinimizerFragment::countDoubleBonds() const
{
    int n = 0;
    for (auto m_bond : m_bonds) {
        if (m_bond->bondOrder == 2) {
            ++n;
        }
    }
    return n;
}
unsigned int sketcherMinimizerFragment::countHeavyAtoms() const
{
    int n = 0;
    for (auto m_atom : m_atoms) {
        if (m_atom->atomicNumber != 6) {
            ++n;
        }
    }
    return n;
}

unsigned int sketcherMinimizerFragment::countConstrainedAtoms() const
{
    int n = 0;
    for (auto atom : m_atoms) {
        if (atom->constrained) {
            ++n;
        }
    }
    return n;
}

unsigned int sketcherMinimizerFragment::countFixedAtoms() const
{
    int n = 0;
    for (auto atom : m_atoms) {
        if (atom->fixed) {
            ++n;
        }
    }
    return n;
}

void sketcherMinimizerFragment::addAtom(sketcherMinimizerAtom* atom)
{
    m_atoms.push_back(atom);
    atom->setFragment(this);
}
void sketcherMinimizerFragment::addBond(sketcherMinimizerBond* bond)
{
    m_bonds.push_back(bond);
}
void sketcherMinimizerFragment::addRing(sketcherMinimizerRing* ring)
{
    m_rings.push_back(ring);
}

void sketcherMinimizerFragment::setAllCoordinatesToTemplate()
{
    for (sketcherMinimizerAtom* atom : m_atoms) {
        atom->setCoordinatesToTemplate();
    }
    if (_bondToParent) {
        _bondToParent->startAtom->setCoordinatesToTemplate();
        _bondToParent->endAtom->setCoordinatesToTemplate();
    }
    for (sketcherMinimizerFragment* child : _children) {
        child->_bondToParent->startAtom->setCoordinatesToTemplate();
        child->_bondToParent->endAtom->setCoordinatesToTemplate();
    }
}

void sketcherMinimizerFragment::storeCoordinateInformation()
{

    _coordinates.clear();
    sketcherMinimizerPointF origin(0.f, 0.f);
    float angle = 0.f;
    if (_bondToParent) {
        origin = _bondToParent->endAtom->getCoordinates();
        angle = atan2(_bondToParent->startAtom->coordinates.y() - origin.y(),
                      -_bondToParent->startAtom->coordinates.x() + origin.x());
    } else {
        if (!constrained && !fixed) {
            origin = m_atoms[0]->getCoordinates();
        }
    }

    float cosine = cos(-angle);
    float sine = sin(-angle);

    for (auto at : m_atoms) {
        sketcherMinimizerPointF c = at->coordinates - origin;
        c.rotate(sine, cosine);
        _coordinates[at] = c;
    }

    for (auto child : _children) {
        sketcherMinimizerAtom* at = child->_bondToParent->endAtom;
        sketcherMinimizerPointF c = at->coordinates - origin;
        c.rotate(sine, cosine);
        _coordinates[at] = c;
    }
}

void sketcherMinimizerFragment::setCoordinates(
    const sketcherMinimizerPointF& position, float angle)
{
    float sine = sin(angle);
    float cosine = cos(angle);
    assert(_coordinates.size() == m_atoms.size() + _children.size());
    for (auto& atom : _coordinates) {
        atom.first->setCoordinates(atom.second);
    }
    for (auto dof : m_dofs) {
        dof->apply();
    }
    for (auto& coords : _coordinates) {
        sketcherMinimizerAtom* atom = coords.first;
        sketcherMinimizerPointF initialCoordinates = atom->getCoordinates();
        initialCoordinates.rotate(sine, cosine);
        atom->setCoordinates(initialCoordinates + position);
    }
}