group.h

#pragma once
#include "core.h"
#include "particle.h"
#include "molecule.h"
#include "units.h"

namespace Faunus {

template <typename T> void inline swap_to_back(T first, T last, T end)
{
    while (end-- > last) {
        std::iter_swap(first++, end);
    }
} 

template <class T> struct IterRange : std::pair<T, T>
{
    using std::pair<T, T>::pair;

    T& begin() { return this->first; }

    T& end() { return this->second; }

    [[nodiscard]] const T& begin() const { return this->first; }

    [[nodiscard]] const T& end() const { return this->second; }

    [[nodiscard]] size_t size() const { return std::distance(this->first, this->second); }

    void resize(size_t n)
    {
        end() += n - size();
        assert(size() == n);
    }

    [[nodiscard]] bool empty() const { return this->first == this->second; }

    void clear()
    {
        this->second = this->first;
        assert(empty());
    }

    [[nodiscard]] std::pair<int, int> to_index(T reference) const
    {
        return {std::distance(reference, begin()), std::distance(reference, end() - 1)};
    } 
}; 

template <class T> class ElasticRange : public IterRange<typename std::vector<T>::iterator>
{
  public:
    using Titer = typename std::vector<T>::iterator;
    using const_iterator = typename std::vector<T>::const_iterator;

  private:
    Titer _trueend;

  public:
    using base = IterRange<typename std::vector<T>::iterator>;
    using base::begin;
    using base::end;
    using base::size;

    ElasticRange(Titer begin, Titer end);
    virtual ~ElasticRange() = default;
    [[nodiscard]] size_t capacity() const;
    [[nodiscard]] auto inactive() const; 
    void
    deactivate(Titer first,
               Titer last); 
    void activate(Titer first, Titer last); 
    Titer& trueend();
    [[nodiscard]] const Titer& trueend() const;
    void relocate(const_iterator oldorigin,
                  Titer neworigin); 
    [[nodiscard]] [[nodiscard]] virtual bool isFull() const;
    [[maybe_unused]] [[nodiscard]] auto numInactive() const; 


    [[maybe_unused]] [[nodiscard]] inline bool
    resizeIsPossible(int number_to_insert_or_delete) const
    {
        auto new_size = static_cast<int>(size()) + number_to_insert_or_delete;
        return (new_size >= 0 && new_size <= capacity());
    }

    inline auto all()
    {
        return std::ranges::subrange(begin(), trueend());
    } 

    [[nodiscard]] inline auto all() const
    {
        return std::ranges::subrange(begin(), trueend());
    } 
};

template <class T> bool ElasticRange<T>::isFull() const
{
    return end() == trueend();
}

template <class T>
ElasticRange<T>::ElasticRange(ElasticRange::Titer begin, ElasticRange::Titer end)
    : base({begin, end})
    , _trueend(end)
{
}

template <class T> size_t ElasticRange<T>::capacity() const
{
    return std::distance(begin(), _trueend);
}

template <class T> auto ElasticRange<T>::inactive() const
{
    return base({end(), _trueend});
}

template <class T>
void ElasticRange<T>::deactivate(ElasticRange::Titer first, ElasticRange::Titer last)
{
    size_t n = std::distance(first, last);
    assert(first >= begin() && last <= end());
    std::rotate(begin(), last, end());
    end() -= n;
    assert(size() + inactive().size() == capacity());
}

template <class T>
void ElasticRange<T>::activate(ElasticRange::Titer first, ElasticRange::Titer last)
{
    size_t n = std::distance(first, last);
    std::rotate(end(), first, _trueend);
    end() += n;
    assert(size() + inactive().size() == capacity());
}

template <class T> typename ElasticRange<T>::Titer& ElasticRange<T>::trueend()
{
    return _trueend;
}

template <class T> const typename ElasticRange<T>::Titer& ElasticRange<T>::trueend() const
{
    return _trueend;
}

template <class T>
void ElasticRange<T>::relocate(ElasticRange::const_iterator oldorigin,
                               ElasticRange::Titer neworigin)
{
    begin() = neworigin + std::distance(oldorigin, const_iterator(begin()));
    end() = neworigin + std::distance(oldorigin, const_iterator(end()));
    trueend() = neworigin + std::distance(oldorigin, const_iterator(trueend()));
}

template <class T> [[maybe_unused]] auto ElasticRange<T>::numInactive() const
{
    return inactive().size();
}

class Group : public ElasticRange<Particle>
{
  public:
    ~Group() override = default;
    using base = ElasticRange<Particle>;
    using iter = typename base::Titer;
    int id = -1;                         
    int conformation_id = 0;             
    Point mass_center = {0.0, 0.0, 0.0}; 

    [[nodiscard]] inline bool isAtomic() const
    {
        return traits().atomic;
    } 

    [[nodiscard]] inline bool isMolecular() const
    {
        return !traits().atomic;
    } 

    [[nodiscard]] bool isFull() const override; 

    std::optional<std::reference_wrapper<Point>>
    massCenter(); 
    [[nodiscard]] std::optional<std::reference_wrapper<const Point>>
    massCenter() const; 

    enum Selectors : unsigned int
    {
        ANY = (1U << 1U),       
        ACTIVE = (1U << 2U),    
        INACTIVE = (1U << 3U),  
        NEUTRAL = (1U << 4U),   
        ATOMIC = (1U << 5U),    
        MOLECULAR = (1U << 6U), 
        FULL = (1U << 7U)       
    };

#pragma clang diagnostic push
#pragma clang diagnostic ignored "-Wc++11-narrowing"

    template <unsigned int mask> [[nodiscard]] bool match() const
    {
        static_assert(mask >= ANY && mask <= FULL);
        if constexpr (mask & ANY) {
            static_assert(mask == ANY, "don't mix ANY with other flags");
            return true;
        }
        if constexpr (mask & ACTIVE) {
            static_assert(!(mask & INACTIVE), "don't mix ACTIVE and INACTIVE");
            if (empty()) {
                return false;
            }
        }
        else if constexpr (mask & INACTIVE) {
            if (!empty()) {
                return false;
            }
        }
        if constexpr (mask & FULL) {
            if (end() != trueend()) {
                return false;
            }
        }
        if constexpr (mask & ATOMIC) {
            static_assert(!(mask & MOLECULAR), "don't mix ATOMIC and MOLECULAR");
            if (isMolecular()) {
                return false;
            }
        }
        else if constexpr (mask & MOLECULAR) {
            if (isAtomic()) {
                return false;
            }
        }
        if constexpr (mask & NEUTRAL) {
            auto _end = (mask & INACTIVE) ? trueend() : end();
            auto _charge = std::accumulate(begin(), _end, 0.0, [](auto sum, const auto& particle) {
                return sum + particle.charge;
            });
            if (std::fabs(_charge) > pc::epsilon_dbl) {
                return false;
            }
        }
        return true;
    }

#pragma clang diagnostic pop

    [[nodiscard]] inline const MoleculeData& traits() const
    {
        assert(id >= 0 && id < Faunus::molecules.size());
        return Faunus::molecules[id];
    } 

    Group(Group& other);
    Group(const Group& other);
    Group(MoleculeData::index_type molid, iter begin, iter end); 
    Group& operator=(const Group& other);   
    Group& shallowCopy(const Group& other); 
    [[nodiscard]] bool contains(const Particle& particle,
                                bool include_inactive = false) const; 
    [[nodiscard]] double mass() const;                                

    auto positions()
    {
        return std::ranges::subrange(begin(), end()) |
               std::views::transform([&](Particle& particle) -> Point& { return particle.pos; });
    } 

    [[nodiscard]] auto positions() const
    {
        return std::ranges::subrange(begin(), end()) |
               std::views::transform(
                   [&](const Particle& particle) -> const Point& { return particle.pos; });
    } 

    [[nodiscard]] AtomData::index_type getParticleIndex(const Particle& particle,
                                                        bool include_inactive = false)
        const; 

    [[nodiscard]] auto findAtomID(AtomData::index_type atomid) const
    {
        return *this |
               std::views::filter([atomid](auto& particle) { return (particle.id == atomid); });
    } 

    inline auto& operator[](size_t index) { return *(begin() + index); }

    inline const auto& operator[](size_t index) const { return *(begin() + index); }

    Particle& at(size_t index);
    [[nodiscard]] const Particle& at(size_t index) const;

    template <std::integral Tint = size_t> auto operator[](std::vector<Tint>& indices)
    {
#ifndef NDEBUG
        if (not indices.empty()) {
            assert(*std::max_element(indices.begin(), indices.end()) < size());
        }
#endif
        return indices |
               std::views::transform([this](auto i) -> Particle& { return *(begin() + i); });
    }

    template <typename TdistanceFunc> void unwrap(const TdistanceFunc& vector_distance)
    {
        if (isMolecular()) {
            for (auto& particle : *this) {
                particle.pos = mass_center + vector_distance(particle.pos, mass_center);
            }
        }
    }

    [[maybe_unused]] void
    wrap(Geometry::BoundaryFunction boundary); 

    void translate(
        const Point& displacement, Geometry::BoundaryFunction boundary = [](Point&) {
        }); 

    void rotate(const Eigen::Quaterniond& quaternion,
                Geometry::BoundaryFunction
                    boundary); 

    void updateMassCenter(Geometry::BoundaryFunction boundary,
                          const Point& approximate_mass_center); 

    void updateMassCenter(Geometry::BoundaryFunction boundary); 

}; 

void to_json(json& j, const Group& group);
void from_json(const json& j, Group& group);

template <unsigned int mask> std::function<bool(const Group&)> getGroupFilter()
{
    return [](const Group& group) { return group.match<mask>(); };
}

template <typename T>
concept RequireGroups =
    std::ranges::range<T> && std::is_convertible_v<std::ranges::range_value_t<T>, Group>;

} // namespace Faunus