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

#include "CoordgenFragmentBuilder.h"
#include "CoordgenMinimizer.h"

#include "sketcherMinimizerFragment.h"
#include "sketcherMinimizer.h"
#include <numeric>
#include <algorithm>
using namespace std;

const int bondLength = BONDLENGTH;
const int PERFECTLY_PLANAR_SYSTEM_SCORE = 50;
const int NON_PLANAR_SYSTEM_SCORE = 1000;
const int UNTREATABLE_SYSTEM_PLANARITY_SCORE = 200000;

void CoordgenFragmentBuilder::initializeCoordinates(
    sketcherMinimizerFragment* fragment) const
{
    assert(fragment->getAtoms().size());
    buildFragment(fragment);
    fragment->storeCoordinateInformation();
}

void CoordgenFragmentBuilder::rotateMainFragment(
    sketcherMinimizerFragment* f) const
{
    if (f->fixed)
        return;
    if (f->isTemplated)
        return;
    if (!f->constrained)
        return;

    int counter = 0;
    sketcherMinimizerPointF constrainOldCenter(0.f, 0.f);
    sketcherMinimizerPointF constrainNewCenter(0.f, 0.f);
    vector<sketcherMinimizerAtom*> fragmentAtoms = f->getAtoms();

    foreach (sketcherMinimizerAtom* a, fragmentAtoms) {
        if (a->constrained) {
            constrainOldCenter += a->templateCoordinates;
            constrainNewCenter += a->coordinates;
            counter++;
        }
    }
    foreach (sketcherMinimizerFragment* child, f->_children) {
        sketcherMinimizerAtom* at = child->_bondToParent->endAtom;
        if (at->constrained) {
            constrainOldCenter += at->templateCoordinates;
            constrainNewCenter += at->coordinates;
            counter++;
        }
    }

    if (counter > 0) {
        constrainOldCenter /= counter;
        constrainNewCenter /= counter;
    }

    vector<sketcherMinimizerPointF> v1, v2;
    foreach (sketcherMinimizerAtom* a, fragmentAtoms) {
        if (a->constrained) {
            v2.push_back(a->coordinates - constrainNewCenter);
            v1.push_back(a->templateCoordinates - constrainOldCenter);
        }
    }
    foreach (sketcherMinimizerFragment* child, f->_children) {
        sketcherMinimizerAtom* at = child->_bondToParent->endAtom;
        if (at->constrained) {
            v2.push_back(at->coordinates - constrainNewCenter);
            v1.push_back(at->templateCoordinates - constrainOldCenter);
        }
    }

    float rotMat[4];
    sketcherMinimizer::alignmentMatrix(v1, v2, rotMat);
    foreach (sketcherMinimizerAtom* a, fragmentAtoms) {
        sketcherMinimizerPointF v = a->getCoordinates() - constrainNewCenter;
        v = sketcherMinimizerPointF(v.x() * rotMat[0] + v.y() * rotMat[1],
                                    v.x() * rotMat[2] + v.y() * rotMat[3]);
        a->setCoordinates(v + constrainOldCenter);
    }
    foreach (sketcherMinimizerFragment* child, f->_children) {
        sketcherMinimizerAtom* at = child->_bondToParent->endAtom;
        sketcherMinimizerPointF v = at->getCoordinates() - constrainNewCenter;
        v = sketcherMinimizerPointF(v.x() * rotMat[0] + v.y() * rotMat[1],
                                    v.x() * rotMat[2] + v.y() * rotMat[3]);
        at->setCoordinates(v + constrainOldCenter);
        at->coordinatesSet = false;
    }
}

bool CoordgenFragmentBuilder::findTemplate(
    const vector<sketcherMinimizerRing*> rings) const
{
    vector<sketcherMinimizerAtom*> allAtoms =
        rings[0]->_atoms[0]->fragment->getAtoms();

    bool foundTemplate = false;
    sketcherMinimizer::loadTemplates();

    vector<sketcherMinimizerAtom*> atoms;
    vector<sketcherMinimizerBond*> bonds;
    map<sketcherMinimizerAtom*, bool> isVisited;

    vector<int> oldIndices;
    foreach (sketcherMinimizerRing* r, rings) {
        foreach (sketcherMinimizerAtom* a, r->_atoms) {
            if (!isVisited[a]) {
                isVisited[a] = true;
                oldIndices.push_back(a->_generalUseN);
                a->_generalUseN = atoms.size();

                atoms.push_back(a);
            }
        }
    }
    foreach (sketcherMinimizerRing* r, rings) {
        foreach (sketcherMinimizerBond* b, r->_bonds) {
            bool found = false;
            foreach (sketcherMinimizerBond* b2, bonds) {
                if (b2 == b) {
                    found = true;
                    break;
                }
            }
            if (!found) {
                bonds.push_back(b);
            }
        }
    }
    vector<unsigned int> mapping;
    for (auto temp : sketcherMinimizer::m_templates.getTemplates()) {
        foundTemplate = sketcherMinimizer::compare(atoms, bonds, temp, mapping);
        if (foundTemplate) {
            if (atoms.size() > 0)
                atoms[0]->fragment->isTemplated = true;
            for (unsigned int i = 0; i < atoms.size(); i++) {
                atoms[i]->setCoordinates(temp->_atoms[mapping[i]]->coordinates *
                                         bondLength);
                atoms[i]->rigid = true;
            }
            foreach (sketcherMinimizerRing* r, rings) {
                r->coordinatesGenerated = true;
            }
            break;
        }
    }
    for (unsigned int i = 0; i < oldIndices.size(); i++) {
        atoms[i]->_generalUseN = oldIndices[i];
    }
    return foundTemplate;
}

sketcherMinimizerRing* CoordgenFragmentBuilder::findCentralRingOfSystem(
    const vector<sketcherMinimizerRing*> rings) const
{
    sketcherMinimizerRing* highest = rings.at(0);
    int high_score = -1;
    foreach (sketcherMinimizerRing* r, rings) {
        int priority = 0;
        if (r->isMacrocycle())
            priority += 1000;
        if (r->_atoms.size() == 6)
            priority += 100;
        priority += r->_atoms.size();
        priority += 10 * (r->fusedWith.size());
        if (priority > high_score) {
            highest = r;
            high_score = priority;
        }
    }
    return highest;
}

void CoordgenFragmentBuilder::generateCoordinatesCentralRings(
    vector<sketcherMinimizerRing*> rings) const
{
    if (rings.size() == 1) { // linearfusedRings
        buildRing(rings[0]);
    } else {
        bool foundTemplate = findTemplate(rings);
        if (!foundTemplate) {
            float planarityScore = newScorePlanarity(rings);
            if (planarityScore < NON_PLANAR_SYSTEM_SCORE) {
                bool needsTemplate = planarityScore > PERFECTLY_PLANAR_SYSTEM_SCORE;
                if (needsTemplate) {
                    findTemplate(rings);
                }
                while (rings.size() != 0) {
                    auto centralRing = findCentralRingOfSystem(rings);
                    buildRing(centralRing);
                    rings.erase(
                        std::remove(rings.begin(), rings.end(), centralRing),
                        rings.end());
                }
                CoordgenMinimizer::maybeMinimizeRings(rings);
            }
            else if (planarityScore > UNTREATABLE_SYSTEM_PLANARITY_SCORE) {
                return;
            }
            else {
                sketcherMinimizerRing* firstRing =
                    findCentralRingOfSystem(rings);
                m_macrocycleBuilder.openCycleAndGenerateCoords(firstRing);
                firstRing->getAtoms()
                    .at(0)
                    ->getMolecule()
                    ->requireMinimization();
            }
        }
    }
}

float CoordgenFragmentBuilder::newScorePlanarity(
    const vector<sketcherMinimizerRing*> rings)
    const // if score > 1000 then it is not planar
{
    float score = 0.f;

    for (auto ring : rings) {
        if (ring->isMacrocycle() &&
            m_macrocycleBuilder.findBondToOpen(ring) == NULL) {
            continue;
        }
        if (ring->isMacrocycle()) {
            for (auto otherRing : ring->fusedWith) {
                if (otherRing->isMacrocycle()) {
                    score += NON_PLANAR_SYSTEM_SCORE;
                }
            }
        }
        for (auto bond : ring->_bonds) {
            if (bond->rings.size() > 2) {
                score += NON_PLANAR_SYSTEM_SCORE * (bond->rings.size() - 2);
            }
        }
        for (auto atom : ring->getAtoms()) {
            if (atom->neighbors.size() > 3) {
                float angle = 0;
                foreach (sketcherMinimizerRing* r, atom->rings) {
                    angle += M_PI - (2 * M_PI / r->_atoms.size());
                }
                if (angle >= 1.99 * M_PI) {
                    score += NON_PLANAR_SYSTEM_SCORE;
                }
            }
        }
    }
    return score;
}

void CoordgenFragmentBuilder::generateCoordinatesSideRings(
    stack<sketcherMinimizerRing*> sideRings) const
{
    while (sideRings.size()) {
        sketcherMinimizerRing* ring = sideRings.top();
        sideRings.pop();
        buildRing(ring);
    }
}

sketcherMinimizerRing*
CoordgenFragmentBuilder::getSharedAtomsWithAlreadyDrawnRing(
    const sketcherMinimizerRing* ring,
    vector<sketcherMinimizerAtom*>& fusionAtoms,
    sketcherMinimizerBond*& fusionBond) const
{
    sketcherMinimizerRing* parent = NULL;
    for (unsigned int i = 0; i < ring->fusedWith.size(); i++) {

        if (ring->fusedWith[i]->coordinatesGenerated) {
            if (!parent)
                parent = ring->fusedWith[i];
            else {
                if (ring->fusedWith[i]->_atoms.size() > parent->_atoms.size())
                    parent = ring->fusedWith[i];
            }
        }
    }
    if (parent) {
        for (unsigned int i = 0; i < parent->fusedWith.size(); i++) {
            if (parent->fusedWith[i] == ring) {
                fusionAtoms = parent->fusionAtoms[i];
            }
        }
        foreach (sketcherMinimizerBond* b, parent->fusionBonds) {
            if (ring->containsAtom(b->startAtom) ||
                ring->containsAtom(b->endAtom)) {
                fusionBond = b;
                break;
            }
        }
    }

    return parent;
}

vector<sketcherMinimizerAtom*> CoordgenFragmentBuilder::orderChainOfAtoms(
    const vector<sketcherMinimizerAtom*> atoms,
    sketcherMinimizerAtom* startAtom)
{
    vector<sketcherMinimizerAtom*> orderedAtoms;
    map<sketcherMinimizerAtom*, bool> stillToAdd;
    foreach (sketcherMinimizerAtom* a, atoms) {
        stillToAdd[a] = true;
    }
    sketcherMinimizerAtom* atomToAdd = startAtom;
    while (orderedAtoms.size() < atoms.size()) {
        orderedAtoms.push_back(atomToAdd);
        stillToAdd[atomToAdd] = false;
        if (orderedAtoms.size() >= atoms.size())
            break;

        foreach (sketcherMinimizerAtom* neighbor, atomToAdd->neighbors) {
            if (stillToAdd[neighbor]) {
                atomToAdd = neighbor;
                break;
            }
        }
    }
    return orderedAtoms;
}

vector<sketcherMinimizerAtom*>
CoordgenFragmentBuilder::orderRingAtoms(const sketcherMinimizerRing* ring)
{
    vector<sketcherMinimizerAtom*> ringAtoms = ring->getAtoms();
    assert(!ringAtoms.empty());
    return orderChainOfAtoms(ringAtoms, ringAtoms.at(0));
}



void CoordgenFragmentBuilder::buildRing(sketcherMinimizerRing* ring) const
{
    if (ring->coordinatesGenerated)
        return;

    vector<sketcherMinimizerAtom*> atoms;

    sketcherMinimizerRing* parent = NULL;
    sketcherMinimizerBond* fusionBond = NULL;
    vector<sketcherMinimizerAtom*> fusionAtoms;
    parent = getSharedAtomsWithAlreadyDrawnRing(ring, fusionAtoms, fusionBond);
    atoms = orderRingAtoms(ring);

    vector<sketcherMinimizerPointF> coords;
    if (atoms.size() >= MACROCYCLE) {
        coords = m_macrocycleBuilder.newMacrocycle(ring, atoms);
    } else {
        coords = listOfCoordinatesFromListofRingAtoms(atoms);
    }
    if (parent) {
        if (fusionAtoms.size() < 2) {
            sketcherMinimizerPointF coordinatesOfPivotAtom;

            sketcherMinimizerAtom* pivotAtom = NULL;
            sketcherMinimizerAtom* pivotAtomOnParent = NULL;
            if (!fusionAtoms.empty()) {
                pivotAtom = fusionAtoms.at(0);
                pivotAtomOnParent = fusionAtoms.at(0);
            } else {
                assert(fusionBond != NULL);
                if (find(atoms.begin(), atoms.end(),
                         fusionBond->getStartAtom()) != atoms.end()) {
                    pivotAtom = fusionBond->getStartAtom();
                    pivotAtomOnParent = fusionBond->getEndAtom();
                } else {
                    pivotAtom = fusionBond->getEndAtom();
                    pivotAtomOnParent = fusionBond->getStartAtom();
                }
            }
            assert(find(atoms.begin(), atoms.end(), pivotAtom) != atoms.end());
            assert(parent->containsAtom(pivotAtomOnParent));

            for (unsigned int i = 0; i < atoms.size(); i++) {
                if (atoms[i] == pivotAtom) {
                    coordinatesOfPivotAtom = coords[i];
                    break;
                }
            }

            for (unsigned int i = 0; i < coords.size(); i++) {
                coords[i] -= coordinatesOfPivotAtom;
            }
            sketcherMinimizerPointF center =
                accumulate(coords.begin(), coords.end(),
                           sketcherMinimizerPointF(0, 0)) /
                coords.size();
            center.normalize();

            int neighborsN = 0;
            sketcherMinimizerPointF neighborsMean(0.f, 0.f);
            foreach (sketcherMinimizerAtom* neighbor,
                     pivotAtomOnParent->neighbors) {
                if (parent->containsAtom(neighbor)) {
                    neighborsMean += neighbor->getCoordinates();
                    ++neighborsN;
                }
            }
            assert(neighborsN == 2);
            neighborsMean *= 0.5;
            sketcherMinimizerPointF pivotAtomCoordinates =
                pivotAtomOnParent->getCoordinates();
            sketcherMinimizerPointF direction =
                pivotAtomCoordinates - neighborsMean;
            direction.normalize();
            sketcherMinimizerPointF bondVector(0.f, 0.f);
            if (fusionBond != NULL) {
                bondVector = direction * bondLength;
            }
            for (unsigned int i = 0; i < coords.size(); i++) {
                sketcherMinimizerPointF coordinates = coords[i];
                coordinates.rotate(center.y(), center.x());
                coordinates.rotate(-direction.y(), direction.x());
                coords[i] = coordinates + pivotAtomCoordinates + bondVector;
            }
            for (unsigned int i = 0; i < coords.size(); i++) {
                atoms[i]->setCoordinates(coords[i]);
            }
            if (fusionBond) {
                // check ZE inversions
                if (!fusionBond->checkStereoChemistry()) {
                    foreach (sketcherMinimizerAtom* atom, atoms) {
                        sketcherMinimizerAtom::mirrorCoordinates(atom,
                                                                 fusionBond);
                    }
                }
            }

        }

        else {
            // build 2 lists of coordinates, mirror of each other.

            vector<sketcherMinimizerPointF> coords2 = coords;
            for (unsigned int i = 0; i < coords2.size(); i++) {
                coords2[i] =
                    sketcherMinimizerPointF(coords2[i].x(), -coords2[i].y());
            }
            map<sketcherMinimizerAtom *, sketcherMinimizerPointF> map1, map2;
            for (unsigned int i = 0; i < atoms.size(); i++) {
                map1[atoms[i]] = coords[i];
                map2[atoms[i]] = coords2[i];
            }
            int lastI = fusionAtoms.size() - 1;
            sketcherMinimizerPointF mean = (fusionAtoms[0]->coordinates +
                                            fusionAtoms[lastI]->coordinates) *
                                           0.5;
            sketcherMinimizerPointF newMean =
                (map1[fusionAtoms[0]] + map1[fusionAtoms[lastI]]) * 0.5;
            sketcherMinimizerPointF v =
                fusionAtoms[lastI]->coordinates - fusionAtoms[0]->coordinates;
            float distance = v.length();
            if (distance < SKETCHER_EPSILON) {
                distance = SKETCHER_EPSILON;
            }
            float sine = v.y() / distance;
            float cosine = v.x() / distance;
            sketcherMinimizerPointF newV =
                (map1[fusionAtoms[lastI]] - map1[fusionAtoms[0]]);
            float newDistance = newV.length();
            float newSine = newV.y() / newDistance;
            float newCosine = newV.x() / newDistance;

            for (unsigned int i = 0; i < atoms.size();
                 i++) { // rotate and translate coords
                sketcherMinimizerPointF d = coords[i] - newMean;
                d.rotate(newSine, newCosine);
                d.rotate(-sine, cosine);
                coords[i] = mean + d;
            }
            newV = (map2[fusionAtoms[lastI]] - map2[fusionAtoms[0]]);
            newMean = (map2[fusionAtoms[0]] + map2[fusionAtoms[lastI]]) * 0.5;

            newDistance = newV.length();
            newSine = newV.y() / newDistance;
            newCosine = newV.x() / newDistance;

            for (unsigned int i = 0; i < atoms.size();
                 i++) { // rotate and translate coords2
                sketcherMinimizerPointF d = coords2[i] - newMean;
                d.rotate(newSine, newCosine);
                d.rotate(-sine, cosine);
                coords2[i] = mean + d;
            }
            // pick the list of coordinates that will place the ring outside the
            // parent i.e further away from its center.
            sketcherMinimizerPointF center = parent->findCenter();
            map<sketcherMinimizerAtom*, bool> fusionMap;
            for (unsigned int i = 0; i < fusionAtoms.size(); i++)
                fusionMap[fusionAtoms[i]] = true;
            float distance1 = 0.f, distance2 = 0.f;
            for (unsigned int i = 0; i < atoms.size(); i++) {
                if (fusionMap[atoms[i]])
                    continue;
                sketcherMinimizerPointF v1 = coords[i] - center;
                sketcherMinimizerPointF v2 = coords2[i] - center;
                distance1 += v1.length();
                distance2 += v2.length();
            }
            vector<sketcherMinimizerPointF>* targetCoords = &coords;
            if (distance1 < distance2)
                targetCoords = &coords2;
            for (unsigned int i = 0; i < atoms.size(); i++) {
                atoms[i]->setCoordinates((*targetCoords)[i]);
            }
            // TODO forbid if double bond to ring
            if (ring->fusedWith.size() == 1 &&
                (*fusionAtoms.begin())->hasNoStereoActiveBonds() &&
                (*fusionAtoms.rbegin())->hasNoStereoActiveBonds()) {
                CoordgenFlipRingDOF* dof =
                    new CoordgenFlipRingDOF(ring, fusionAtoms);
                ring->getAtoms().at(0)->fragment->addDof(dof);
            }
        }

    } else {
        for (unsigned int i = 0; i < coords.size(); i++) {
            atoms[i]->setCoordinates(coords[i]);
        }
    }
    ring->coordinatesGenerated = true;
}

vector<sketcherMinimizerPointF>
CoordgenFragmentBuilder::listOfCoordinatesFromListofRingAtoms(
    const vector<sketcherMinimizerAtom*> atoms)
{
    vector<sketcherMinimizerPointF> out;
    assert(atoms.size());
    float a = 2 * M_PI / atoms.size();
    sketcherMinimizerPointF coords(0.f, 0.f);
    float angle = 0;
    for (unsigned int n = 0; n < atoms.size(); n++) {
        out.push_back(coords);
        angle = a * n;
        coords += sketcherMinimizerPointF(cos(angle) * bondLength,
                                          -sin(angle) * bondLength);
    }
    return out;
}

void CoordgenFragmentBuilder::buildNonRingAtoms(
    sketcherMinimizerFragment* fragment) const
{
    set<sketcherMinimizerAtom*> isAtomVisited;
    queue<sketcherMinimizerAtom*> q;
    for (sketcherMinimizerAtom* atom : fragment->getAtoms()) {
        if (!atom->getRings().empty()) {
            q.push(atom);
            isAtomVisited.insert(atom);
        }
    }
    if (q.empty()) {
        if (fragment->getParent()) {
            fragment->_bondToParent->getStartAtom()->setCoordinates(
                sketcherMinimizerPointF(-bondLength, 0));
        }
        sketcherMinimizerAtom* atomToStartFrom =
            (fragment->getParent() ? fragment->_bondToParent->getEndAtom()
                                   : fragment->getAtoms()[0]);
        atomToStartFrom->setCoordinates(sketcherMinimizerPointF(0, 0));
        q.push(atomToStartFrom);
        isAtomVisited.insert(atomToStartFrom);
    }

    while (q.size()) {
        generateCoordinatesNeighborsOfFirstAtomInQueue(q, isAtomVisited,
                                                       fragment);
    }
}

vector<float> CoordgenFragmentBuilder::neighborsAnglesAtCenter(
    const sketcherMinimizerAtom* atom) const
{

    int angleDivision = atom->neighbors.size();
    vector<float> angles;
    if (!m_evenAngles) {
        if (atom->neighbors.size() == 2) {
            if (atom->getAtomicNumber() == 6 ||
                !atom->neighbors.at(0)->crossLayout ||
                !atom->neighbors.at(1)->crossLayout) {
                angleDivision = 3;
            }
            if (atom->bonds[0]->getBondOrder() +
                    atom->bonds[1]->getBondOrder() >=
                4)
                angleDivision = 2;

        } else if (atom->neighbors.size() == 4 && !atom->crossLayout) {
            angles.push_back(M_PI / 3);
            angles.push_back(M_PI * 0.5);
            angles.push_back(2 * M_PI / 3);
            angles.push_back(M_PI * 0.5);
        }
    }
    if (!angles.size()) {
        for (unsigned int i = 0; i < atom->neighbors.size(); i++) {
            angles.push_back(2 * M_PI / angleDivision);
        }
    }
    return angles;
}

void CoordgenFragmentBuilder::initializeVariablesForNeighboursCoordinates(
    sketcherMinimizerAtom* atom, set<sketcherMinimizerAtom*>& isAtomVisited,
    sketcherMinimizerPointF& startCoordinates,
    vector<sketcherMinimizerAtom*>& orderedNeighbors,
    vector<float>& angles) const
{
    if (atom->rings.size()) {
        return initializeVariablesForNeighboursCoordinatesRingAtom(
            atom, isAtomVisited, startCoordinates, orderedNeighbors, angles);
    }
    orderedNeighbors.clear();
    if (atom->neighbors.size() != 4) {
        orderedNeighbors = atom->neighbors;
    } else {
        vector<sketcherMinimizerAtomPriority> atomPriorities;
        for (unsigned int i = 0; i < atom->neighbors.size(); i++) {
            sketcherMinimizerAtomPriority p;
            p.a = atom->neighbors[i];
            atomPriorities.push_back(p);
        }
        sketcherMinimizerAtom::orderAtomPriorities(atomPriorities, atom);
        foreach (sketcherMinimizerAtomPriority ap, atomPriorities) {
            orderedNeighbors.push_back(ap.a);
        }
    }
    unsigned int startN = 0;
    for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
        if (isAtomVisited.find(orderedNeighbors[i]) != isAtomVisited.end()) {
            startCoordinates =
                orderedNeighbors[i]->getCoordinates() - atom->getCoordinates();
            startN = i;
            break;
        }
    }
    for (unsigned int i = 0; i < startN; i++) {
        orderedNeighbors.push_back(orderedNeighbors.at(0));
        orderedNeighbors.erase(orderedNeighbors.begin());
    }
    angles = neighborsAnglesAtCenter(atom);
}

void CoordgenFragmentBuilder::
    initializeVariablesForNeighboursCoordinatesRingAtom(
        const sketcherMinimizerAtom* atom,
        set<sketcherMinimizerAtom*>& ,
        sketcherMinimizerPointF& startCoordinates,
        vector<sketcherMinimizerAtom*>& orderedNeighbors,
        vector<float>& angles) const
{

    vector<pair<float, sketcherMinimizerAtom*>> ringNeighboursAndAngles;
    orderedNeighbors.clear();
    foreach (sketcherMinimizerAtom* neigh, atom->neighbors) {
        if (sketcherMinimizer::sameRing(neigh, atom)) {
            float ang = atan2(neigh->coordinates.y() - atom->coordinates.y(),
                              neigh->coordinates.x() - atom->coordinates.x());
            if (ang < 0)
                ang += 2 * M_PI;
            pair<float, sketcherMinimizerAtom*> pairToAdd(ang, neigh);
            ringNeighboursAndAngles.push_back(pairToAdd);
        } else {
            orderedNeighbors.push_back(neigh);
        }
    }
    stable_sort(ringNeighboursAndAngles.begin(), ringNeighboursAndAngles.end());
    vector<float> gaps;
    vector<float> scaledGaps;
    for (unsigned int i = 0; i < ringNeighboursAndAngles.size(); i++) {
        int next = (i + 1) % ringNeighboursAndAngles.size();
        float gap = ringNeighboursAndAngles[next].first -
                    ringNeighboursAndAngles[i].first;
        if (gap < 0)
            gap += 2 * M_PI;
        gaps.push_back(gap);
        bool rin = false;
        float mid = ringNeighboursAndAngles[i].first + gap * 0.5;
        float sine = sin(-mid);
        float cosine = cos(-mid);
        sketcherMinimizerPointF p(bondLength * 0.1, 0);
        p.rotate(sine, cosine);
        sketcherMinimizerPointF point = atom->coordinates + p;
        vector<sketcherMinimizerRing*> rings = atom->getFragment()->getRings();
        foreach (sketcherMinimizerRing* r, rings) {
            if (r->isMacrocycle())
                continue;
            if (r->contains(point)) {
                rin = true;
                break;
            }
        }
        float scaledGap = gap;
        if (gap > M_PI)
            scaledGap *= 10;
        else if (rin)
            scaledGap *= 0.2;
        scaledGaps.push_back(scaledGap);
    }
    int bestI = 0;
    for (unsigned int i = 0; i < gaps.size(); i++) {
        if (scaledGaps[i] > scaledGaps[bestI]) {
            bestI = i;
        }
    }
    angles.clear();
    float biggestGap = gaps[bestI];
    startCoordinates = ringNeighboursAndAngles[bestI].second->getCoordinates() -
                       atom->getCoordinates();
    for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
        angles.push_back(-biggestGap / (orderedNeighbors.size() + 1));
    }
}

void CoordgenFragmentBuilder::generateCoordinatesNeighborsOfFirstAtomInQueue(
    queue<sketcherMinimizerAtom*>& atomQueue,
    set<sketcherMinimizerAtom*>& isAtomVisited,
    const sketcherMinimizerFragment* fragment) const
{

    sketcherMinimizerAtom* at = atomQueue.front();
    atomQueue.pop();
    sketcherMinimizerPointF initialCoords(-bondLength, 0.f);
    vector<sketcherMinimizerAtom*> orderedNeighbors;
    vector<float> angles;

    initializeVariablesForNeighboursCoordinates(
        at, isAtomVisited, initialCoords, orderedNeighbors, angles);

    for (unsigned int i = 0; i < orderedNeighbors.size(); i++) {
        int n = i;
        sketcherMinimizerAtom* neigh = orderedNeighbors[n];
        if (isAtomVisited.find(neigh) != isAtomVisited.end())
            continue;
        float s = sin(angles[n]);
        float c = cos(angles[n]);
        initialCoords.rotate(s, c);
        neigh->setCoordinates(at->coordinates + initialCoords);
        if (neigh->fragment == fragment) {
            isAtomVisited.insert(neigh);
            atomQueue.push(neigh);
        } else {
            neigh->coordinatesSet = false;
        }
        if (at->needsCheckForClashes) {
            neigh->needsCheckForClashes = true;
            CoordgenMinimizer::checkForClashes(neigh);
        }
        for (auto dof : at->fragment->getDofsOfAtom(at)) {
            if (dof->getFragment() == fragment) {
                dof->addAtom(neigh);
            }
        }
    }
    avoidZEInversions(at, isAtomVisited);
    maybeAddMacrocycleDOF(at);
    for (auto neighbor : at->neighbors) {
        if (!sketcherMinimizerAtom::shareARing(at, neighbor) &&
            at->getFragment() == neighbor->getFragment()) {
            CoordgenScaleAtomsDOF* dof = new CoordgenScaleAtomsDOF(at);
            dof->addAtom(neighbor);
            at->fragment->addDof(dof);
        }
    }
}

void CoordgenFragmentBuilder::maybeAddMacrocycleDOF(
    sketcherMinimizerAtom* atom) const
{
    if (atom->getRings().size() == 1 &&
        atom->getRings().at(0)->isMacrocycle() && atom->neighbors.size() == 3) {
        for (auto bond : atom->getBonds()) {
            if (bond->isStereo() && !bond->isTerminal())
                return;
        }
        for (auto neighbor : atom->neighbors) {
            if (!sketcherMinimizerAtom::shareARing(atom, neighbor)) {
                CoordgenInvertBondDOF* dof =
                    new CoordgenInvertBondDOF(atom, neighbor);
                atom->fragment->addDof(dof);
            }
        }
    }
}

void CoordgenFragmentBuilder::avoidZEInversions(
    const sketcherMinimizerAtom* at,
    set<sketcherMinimizerAtom*>& isAtomVisited) const
{
    if (!at->getRings().empty())
        return;
    sketcherMinimizerBond* doubleBond = NULL;
    vector<sketcherMinimizerAtom*> atomsToMirror;
    for (unsigned int i = 0; i < at->bonds.size(); i++) {
        if (at->bonds[i]->isStereo() &&
            isAtomVisited.find(at->neighbors[i]) != isAtomVisited.end()) {
            doubleBond = at->bonds[i];
        } else {
            atomsToMirror.push_back(at->neighbors[i]);
        }
    }
    if (doubleBond == NULL)
        return;
    if (atomsToMirror.size() && doubleBond) {
        sketcherMinimizerAtom* firstCIPNeighborStart =
            doubleBond->startAtomCIPFirstNeighbor();
        if (firstCIPNeighborStart == NULL)
            return;
        sketcherMinimizerAtom* firstCIPNeighborEnd =
            doubleBond->endAtomCIPFirstNeighbor();
        if (firstCIPNeighborEnd == NULL)
            return;
        if (!doubleBond->checkStereoChemistry()) {
            foreach (sketcherMinimizerAtom* a, atomsToMirror) {
                sketcherMinimizerAtom::mirrorCoordinates(a, doubleBond);
            }
        }
    }
}

void CoordgenFragmentBuilder::buildFragment(
    sketcherMinimizerFragment* fragment) const
{
    buildRings(fragment);
    buildNonRingAtoms(fragment);
    CoordgenMinimizer::avoidInternalClashes(fragment);
    fallbackIfNanCoordinates(fragment);
    if (!fragment->getParent() && fragment->constrained)
        rotateMainFragment(fragment);

    // reset coordinates for fixed fragments:
    if (fragment->fixed) {
        fragment->setAllCoordinatesToTemplate();
    }
}


void CoordgenFragmentBuilder::fallbackIfNanCoordinates(
    sketcherMinimizerFragment* fragment) const
{
    vector<sketcherMinimizerAtom*> fragmentAtoms = fragment->getAtoms();
    if (CoordgenMinimizer::hasNaNCoordinates(fragmentAtoms) &&
        CoordgenMinimizer::hasValid3DCoordinates(fragmentAtoms)) {
        CoordgenMinimizer::fallbackOn3DCoordinates(fragmentAtoms);
    }
}

void CoordgenFragmentBuilder::initializeFusedRingInformation(
    sketcherMinimizerFragment* fragment) const
{
    /*
     for each ring in the fragment find out which other rings it shares atoms
     with,
     or if they are connected in other ways (i.e. via a double bond).
     */
    if (fragment->getRings().size() < 2)
        return;
    vector<sketcherMinimizerAtom*> fragmentAtoms = fragment->getAtoms();
    foreach (sketcherMinimizerAtom* atom, fragmentAtoms) {
        if (atom->rings.size() > 1) {
            for (unsigned int i = 0; i < atom->rings.size(); i++) {
                for (unsigned int j = i + 1; j < atom->rings.size(); j++) {

                    bool found = false;
                    for (unsigned int r = 0;
                         r < atom->rings[i]->fusedWith.size(); r++) {
                        if (atom->rings[i]->fusedWith[r] == atom->rings[j]) {
                            atom->rings[i]->fusionAtoms[r].push_back(atom);
                            for (unsigned int rr = 0;
                                 rr < atom->rings[j]->fusedWith.size(); rr++) {

                                if (atom->rings[j]->fusedWith[rr] ==
                                    atom->rings[i]) {
                                    atom->rings[j]->fusionAtoms[rr].push_back(
                                        atom);
                                    break;
                                }
                            }
                            found = true;
                            break;
                        }
                    }
                    if (!found) {
                        atom->rings[i]->fusedWith.push_back(atom->rings[j]);
                        atom->rings[j]->fusedWith.push_back(atom->rings[i]);
                        vector<sketcherMinimizerAtom*> ats;
                        ats.push_back(atom);
                        atom->rings[i]->fusionAtoms.push_back(ats);
                        atom->rings[j]->fusionAtoms.push_back(ats);
                    }
                }
            }
        }
    }
    vector<sketcherMinimizerBond*> bonds = fragment->getBonds();
    foreach (sketcherMinimizerBond* b, bonds) {
        if ((b->bondOrder != 1) && (b->startAtom->rings.size()) &&
            (b->endAtom->rings.size()) &&
            (!sketcherMinimizer::sameRing(b->startAtom, b->endAtom))) {
            foreach (sketcherMinimizerRing* r, b->startAtom->rings) {
                foreach (sketcherMinimizerRing* r2, b->endAtom->rings) {
                    r->fusedWith.push_back(r2);
                    r2->fusedWith.push_back(r);
                    vector<sketcherMinimizerAtom*> ats;
                    r->fusionAtoms.push_back(ats);
                    r2->fusionAtoms.push_back(ats);
                    r->fusionBonds.push_back(b);
                    r2->fusionBonds.push_back(b);
                }
            }
        }
    }
    vector<sketcherMinimizerRing*> fragmentRings = fragment->getRings();
    foreach (sketcherMinimizerRing* r, fragmentRings) {
        for (unsigned int i = 0; i < r->fusedWith.size(); i++) {
            vector<sketcherMinimizerAtom*> fusionAtoms = r->fusionAtoms[i];
            if (fusionAtoms.size() <= 2)
                continue;
            // find an atom with only one neighbor in the vector
            sketcherMinimizerAtom* startAtom = NULL;
            foreach (sketcherMinimizerAtom* a, fusionAtoms) {
                int counter = 0;
                foreach (sketcherMinimizerAtom* n, a->neighbors) {
                    if (find(fusionAtoms.begin(), fusionAtoms.end(), n) !=
                        fusionAtoms.end()) {
                        counter++;
                    }
                }
                if (counter == 1) {
                    startAtom = a;
                    break;
                }
            }
            assert(startAtom != NULL);
            r->fusionAtoms[i] = orderChainOfAtoms(fusionAtoms, startAtom);
        }
    }

    foreach (sketcherMinimizerRing* r, fragmentRings) {
        for (unsigned int ii = 0; ii < r->fusedWith.size(); ii++) {
            if (r->fusionAtoms[ii].size() > 2) {
                foreach (sketcherMinimizerAtom* a,
                         r->fusionAtoms[ii]) { // looking for atoms with no
                                               // neighbors in r but not in
                                               // r->fusedWith[ii]
                    bool found = false;
                    foreach (sketcherMinimizerAtom* n, a->neighbors) {
                        int nn = 0;
                        foreach (sketcherMinimizerRing* nr, n->rings) {
                            if (nr == r)
                                nn++;
                            else if (nr == r->fusedWith[ii])
                                nn++;
                        }
                        if (nn == 1) {
                            found = true;
                            break;
                        }
                    }
                    if (!found) {
                        a->isSharedAndInner = true;
                    }
                }
            }
        }
    }
}

void CoordgenFragmentBuilder::simplifyRingSystem(
    const vector<sketcherMinimizerRing*> allRings,
    stack<sketcherMinimizerRing*>& sideRings,
    vector<sketcherMinimizerRing*>& centralRings) const
{
    map<sketcherMinimizerRing*, bool> alreadyChosenAsSide;

    bool found = true;
    /*simplify the structure by iteratively removing rings that have only one
     fusion partner. These side rings will have coordinates generated
     separately,
     one at a time, starting from the closest to the core.
     If we think of rings as atoms and fusion between two rings as bonds, this
     is the same as stripping away non-ring atoms from a fragment to generate
     coordinate
     for the ring system first
     */
    while (found) {
        found = false;
        foreach (sketcherMinimizerRing* r, allRings) {
            if (alreadyChosenAsSide[r])
                continue;
            if (r->isMacrocycle())
                continue;
            int n = 0;
            for (auto ringCounter = 0u; ringCounter < r->fusedWith.size();
                 ++ringCounter) {
                auto fusedRing = r->fusedWith.at(ringCounter);
                if (!alreadyChosenAsSide[fusedRing]) {
                    n++;
                    if (fusedRing->isMacrocycle()) {
                        /*disable macrocycles*/
                        n++;
                    }
                    if (r->fusionAtoms.at(ringCounter).size() > 3) {
                        /*disable rings that share too many atoms*/
                        n++;
                    }
                }
            }
            if (n == 1) /*ring is fused with only one ring,
                         it will be removed from the core of the system
                         and treated as a side ring
                         */
            {
                sideRings.push(r);
                found = true;
                alreadyChosenAsSide[r] = true;
                r->side = true;
            }
        }
    }
    /*
     everything that is left after removing the side rings is the core of the
     ring system
     */
    foreach (sketcherMinimizerRing* ring, allRings) {
        if (!alreadyChosenAsSide[ring])
            centralRings.push_back(ring);
    }
}

void CoordgenFragmentBuilder::buildRings(
    sketcherMinimizerFragment* fragment) const
{
    if (fragment->getRings().empty())
        return;
    initializeFusedRingInformation(fragment);
    vector<sketcherMinimizerRing*> fragmentRings = fragment->getRings();
    if (fragmentRings.size() == 1) {
        buildRing(fragmentRings[0]);
    } else { // more than 1 ring
        stack<sketcherMinimizerRing*> sideRings;
        vector<sketcherMinimizerRing*> centralRings;
        simplifyRingSystem(fragmentRings, sideRings, centralRings);
        generateCoordinatesCentralRings(centralRings);
        generateCoordinatesSideRings(sideRings);
    }
}