space.h

#pragma once
#include "core.h"
#include "geometry.h"
#include "group.h"
#include "molecule.h"
#include <range/v3/view/join.hpp>
#include <range/v3/view/filter.hpp>
#include <range/v3/range/conversion.hpp>

namespace Faunus {

struct Change
{
    using index_type = std::size_t;
    bool everything = false;    
    bool volume_change = false; 
    bool matter_change = false; 
    bool moved_to_moved_interactions =
        true; 
    bool disable_translational_entropy = false; 

    struct GroupChange
    {
        index_type group_index; 
        bool dNatomic = false;  
        bool dNswap = false;    
        bool internal = false;  
        bool all =
            false; 
        std::vector<index_type>
            relative_atom_indices; 

        bool
        operator<(const GroupChange& other) const; 
        void sort();                               
    };

    std::vector<GroupChange> groups; 

    [[nodiscard]] std::optional<std::pair<index_type, index_type>> singleParticleChange() const;

    [[nodiscard]] inline auto touchedGroupIndex() const
    {
        return std::views::transform(groups, &GroupChange::group_index);
    }

    [[maybe_unused]] [[nodiscard]] std::vector<index_type>
    touchedParticleIndex(const std::vector<Group>&) const;

    void clear();                     
    [[nodiscard]] bool empty() const; 
    explicit operator bool() const;   
    void sanityCheck(
        const std::vector<Group>& group_vector) const; 
} __attribute__((aligned(32)));

void to_json(json& j, const Change::GroupChange& group_change); 
void to_json(json& j, const Change& change);                    

class Space
{
  public:
    using GeometryType = Geometry::Chameleon;
    using GroupType = Group; 
    using GroupVector = std::vector<GroupType>;
    using GroupRefVector = std::vector<std::reference_wrapper<Group>>;
    using ConstGroupRefVector = std::vector<std::reference_wrapper<const Group>>;
    using ScaleVolumeTrigger = std::function<void(Space&, double, double)>;
    using ChangeTrigger = std::function<void(Space&, const Change&)>;
    using SyncTrigger = std::function<void(Space&, const Space&, const Change&)>;

  private:
    std::map<MoleculeData::index_type, std::size_t> implicit_reservoir;

    std::vector<ChangeTrigger> changeTriggers; 
    std::vector<SyncTrigger>
        onSyncTriggers; 

  public:
    ParticleVector particles; 
    GroupVector groups;       
    GeometryType geometry;    
    std::vector<ScaleVolumeTrigger>
        scaleVolumeTriggers; 

    [[nodiscard]] const std::map<MoleculeData::index_type, std::size_t>&
    getImplicitReservoir() const;                                            
    std::map<MoleculeData::index_type, std::size_t>& getImplicitReservoir(); 

    enum class Selection
    {
        ALL,
        ACTIVE,
        INACTIVE,
        ALL_NEUTRAL,
        ACTIVE_NEUTRAL,
        INACTIVE_NEUTRAL
    };

    void clear(); 
    GroupType& addGroup(MoleculeData::index_type molid,
                        const ParticleVector& particles); 
    GroupVector::iterator findGroupContaining(
        const Particle& particle,
        bool include_inactive = false); 
    GroupVector::iterator
    findGroupContaining(AtomData::index_type atom_index); 
    [[nodiscard]] size_t
    numParticles(Selection selection = Selection::ACTIVE) const; 

    Point
    scaleVolume(double new_volume,
                Geometry::VolumeMethod method =
                    Geometry::VolumeMethod::ISOTROPIC); 

    GroupVector::iterator
    randomMolecule(MoleculeData::index_type molid, Random& rand,
                   Selection selection = Selection::ACTIVE); 

    json info();

    [[nodiscard]] std::size_t
    getGroupIndex(const GroupType& group) const; 
    [[nodiscard]] std::size_t
    getFirstParticleIndex(const GroupType& group) const; 
    [[nodiscard]] std::size_t getFirstActiveParticleIndex(
        const GroupType& group) const; 

    template <std::forward_iterator iterator,
              class copy_operation = std::function<void(const Particle&, Particle&)>>
    void updateParticles(
        const iterator begin, const iterator end, ParticleVector::iterator destination,
        copy_operation copy_function = [](const Particle& src, Particle& dst) { dst = src; })
    {

        const auto size = std::distance(begin, end); // number of affected particles

        assert(destination >= particles.begin() && destination < particles.end());
        assert(size <= std::distance(destination, particles.end()));

        auto affected_groups =
            groups | std::views::filter([=](auto& group) {
                return (group.begin() < destination + size) && (group.end() > destination);
            }); // filtered group with affected groups only. Note we copy in org. `destination`

        // copy data from source range (this modifies `destination`)
        std::for_each(begin, end,
                      [&](const auto& source) { copy_function(source, *destination++); });

        std::ranges::for_each(affected_groups, [&](Group& group) {
            group.updateMassCenter(geometry.getBoundaryFunc(), group.begin()->pos);
        });
    }

    void updateInternalState(const Change& change);

    [[nodiscard]] auto positions() const
    {
        return std::views::transform(particles,
                                     [](auto& particle) -> const Point& { return particle.pos; });
    }

    auto positions() { return std::views::transform(particles, &Particle::pos); }

    static std::function<bool(const GroupType&)> getGroupFilter(MoleculeData::index_type molid,
                                                                const Selection& selection)
    {
        auto is_active = [](const GroupType& group) { return group.size() == group.capacity(); };

        auto is_neutral = [](RequireParticleIterator auto begin, RequireParticleIterator auto end) {
            auto charge = std::accumulate(
                begin, end, 0.0, [](auto sum, auto& particle) { return sum + particle.charge; });
            return (fabs(charge) < 1e-6);
        }; 

        std::function<bool(const GroupType&)> f; 

        switch (selection) {
        case (Selection::ALL):
            f = []([[maybe_unused]] auto& group) { return true; };
            break;
        case (Selection::INACTIVE):
            f = [=](auto& group) { return !is_active(group); };
            break;
        case (Selection::ACTIVE):
            f = [=](auto& group) { return is_active(group); };
            break;
        case (Selection::ALL_NEUTRAL):
            f = [=](auto& group) { return is_neutral(group.begin(), group.trueend()); };
            break;
        case (Selection::INACTIVE_NEUTRAL):
            f = [=](auto& group) {
                return !is_active(group) && is_neutral(group.begin(), group.trueend());
            };
            break;
        case (Selection::ACTIVE_NEUTRAL):
            f = [=](auto& group) {
                return is_active(group) && is_neutral(group.begin(), group.end());
            };
            break;
        }
        return [f, molid](auto& group) { return group.id == molid && f(group); };
    }

    auto findMolecules(MoleculeData::index_type molid, Selection selection = Selection::ACTIVE)
    {
        auto group_filter = getGroupFilter(molid, selection);
        return groups | ranges::cpp20::views::filter(group_filter);
    }

    [[nodiscard]] auto findMolecules(MoleculeData::index_type molid,
                                     Selection selection = Selection::ACTIVE) const
    {
        auto group_filter = getGroupFilter(molid, selection);
        return groups | ranges::cpp20::views::filter(group_filter);
    }

    auto activeParticles()
    {
        return groups | ranges::cpp20::views::join;
    } 

    [[nodiscard]] auto activeParticles() const
    {
        return groups | ranges::cpp20::views::join;
    } 

    template <std::integral index_type = int>
    auto toIndices(const RequireParticles auto& particle_range) const
    {
        auto to_index = [&](auto& particle) {
            const auto index = std::addressof(particle) - std::addressof(particles.at(0));
            if (index < 0 || index >= particles.size()) {
                throw std::out_of_range("particle range outside Space");
            }
            return static_cast<index_type>(index);
        };
        return particle_range | std::views::transform(to_index) | ranges::to_vector;
    }

    auto findAtoms(AtomData::index_type atomid)
    {
        return activeParticles() | std::views::filter([atomid](const Particle& particle) {
                   return particle.id == atomid;
               });
    }

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

    [[nodiscard]] size_t countAtoms(AtomData::index_type atomid) const; 

    template <unsigned int mask> auto numMolecules(MoleculeData::index_type molid) const
    {
        auto filter = [&](const GroupType& group) {
            return (group.id == molid) ? group.template match<mask>() : false;
        };
        return std::count_if(groups.begin(), groups.end(), filter);
    }

    void sync(const Space& other,
              const Change& change); 

}; // end of space

void to_json(json& j, const Space& spc);   
void from_json(const json& j, Space& spc); 

class InsertMoleculesInSpace
{
  private:
    static void insertAtomicGroups(MoleculeData& moldata, Space& spc, size_t num_molecules,
                                   size_t num_inactive_molecules);

    static void insertMolecularGroups(MoleculeData& moldata, Space& spc, size_t num_molecules,
                                      size_t num_inactive);

    static void setPositionsForTrailingGroups(Space& spc, size_t num_molecules,
                                              const Faunus::ParticleVector&, const Point&);

    static void insertImplicitGroups(const MoleculeData& moldata, Space& spc, size_t num_molecules);

    static size_t getNumberOfMolecules(const json& j, double volume,
                                       const std::string& molecule_name);

    static size_t getNumberOfInactiveMolecules(const json& j, size_t number_of_molecules);

    static ParticleVector getExternalPositions(const json& j, const std::string& molname);

    static void insertItem(const std::string& molname, const json& properties, Space& spc);

    static void reserveMemory(const json& j, Space& spc);

  public:
    static void insertMolecules(const json& j, Space& spc);
}; // end of insertMolecules class

namespace SpaceFactory {

void makeNaCl(Space& space, size_t num_particles,
              const Geometry::Chameleon& geometry); 
void makeWater(Space& space, size_t num_particles,
               const Geometry::Chameleon& geometry); 

} // namespace SpaceFactory

} // namespace Faunus