analysis.h

#pragma once

#include "space.h"
#include "io.h"
#include "scatter.h"
#include "reactioncoordinate.h"
#include "aux/timers.h"
#include "aux/table_2d.h"
#include "aux/equidistant_table.h"
#include "aux/sparsehistogram.h"
#include <Eigen/SparseCore>
#include <limits>
#include <memory>
#include <set>
#include <string_view>
#include <string>

namespace Faunus::pairpotential {
class NewCoulombGalore;
}

namespace Faunus::Energy {
class Hamiltonian;
class EnergyTerm;
class Penalty;
} // namespace Faunus::Energy

namespace Faunus::SASA {
class SASABase;
}

namespace Faunus::analysis {

class Analysis
{
  private:
    virtual void _to_json(json&) const;   
    virtual void _from_json(const json&); 
    virtual void _sample() = 0;           
    virtual void _to_disk();              
    int number_of_steps = 0;              
    int number_of_skipped_steps = 0;      
    TimeRelativeOfTotal<std::chrono::microseconds> timer; 

  protected:
    const Space& spc;          
    int sample_interval = 0;   
    int number_of_samples = 0; 

  public:
    const std::string name;                     
    std::string cite;                           
    void to_json(json& j) const;                
    void from_json(const json& j);              
    void to_disk();                             
    void sample();                              
    [[nodiscard]] int getNumberOfSteps() const; 
    Analysis(const Space& spc, std::string_view name);
    Analysis(const Space& spc, std::string_view name, int sample_interval,
             int number_of_skipped_steps);
    virtual ~Analysis() = default;
};

void to_json(json& j, const Analysis& base);

std::unique_ptr<Analysis> createAnalysis(const std::string& name, const json& j, Space& spc,
                                         Energy::Hamiltonian& pot);

class CombinedAnalysis : public BasePointerVector<Analysis>
{
  public:
    CombinedAnalysis(const json& json_array, Space& spc, Energy::Hamiltonian& pot);
    void sample();
    void to_disk(); 
};

class PerturbationAnalysis : public Analysis
{
  private:
    void _to_disk() override;

  protected:
    Space& mutable_space; 
    Energy::EnergyTerm& pot;
    std::string filename;                             
    std::unique_ptr<std::ostream> stream = nullptr;   
    Change change;                                    
    Average<double> mean_exponentiated_energy_change; 
    bool collectWidomAverage(double energy_change);   
    PerturbationAnalysis(const std::string& name, Energy::EnergyTerm& pot, Space& spc,
                         const std::string& filename = ""s);
    [[nodiscard]] double
    meanFreeEnergy() const; 
};

class FileReactionCoordinate : public Analysis
{
  private:
    Average<double> mean_reaction_coordinate;
    const std::string filename;
    std::unique_ptr<std::ostream> stream = nullptr;
    const std::unique_ptr<ReactionCoordinate::ReactionCoordinateBase> reaction_coordinate;

    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;
    FileReactionCoordinate(
        const Space& spc, const std::string& filename,
        std::unique_ptr<ReactionCoordinate::ReactionCoordinateBase> reaction_coordinate);

  public:
    FileReactionCoordinate(const json& j, const Space& spc);
};

class AtomicDisplacement : public Analysis
{
  protected:
    MoleculeData::index_type molid;

  private:
    using average_type = Average<double>;
    std::unique_ptr<std::ostream>
        single_position_stream; 
    std::vector<Point> reference_positions;    
    std::vector<Point> previous_positions;     
    std::vector<Eigen::Vector3i> cell_indices; 
    std::vector<average_type>
        mean_squared_displacement;    
    double max_possible_displacement; 

    SparseHistogram<double> displacement_histogram; 
    std::string displacement_histogram_filename;    
    int reference_reset_interval =
        std::numeric_limits<int>::max(); 

    [[nodiscard]] virtual PointVector
    getPositions() const; 
    void resetReferencePosition(const Point& position,
                                int index); 
    Point getOffset(const Point& diff, Eigen::Vector3i& cell) const; 
    void sampleDisplacementFromReference(const Point& position, int index);
    void _sample() override;
    void _to_json(json& j) const override;
    void _to_disk() override;

  public:
    AtomicDisplacement(const json& j, const Space& spc, std::string_view name = "displacement");
};

class MassCenterDisplacement : public AtomicDisplacement
{
  private:
    [[nodiscard]] PointVector
    getPositions() const override; 
  public:
    MassCenterDisplacement(const json& j, const Space& spc,
                           std::string_view name = "displacement_com");
};

class WidomInsertion : public PerturbationAnalysis
{
    std::shared_ptr<MoleculeInserter> inserter; 
    int number_of_insertions;                   
    MoleculeData::index_type molid;             
    bool absolute_z_coords = false;             

    void selectGhostGroup(); 
    void updateGroup(Space::GroupType& group, const ParticleVector& particles);
    void _sample() override; 
    void _to_json(json& j) const override;
    void _from_json(const json& j) override;

  public:
    WidomInsertion(const json& j, Space& spc, Energy::Hamiltonian& pot);
};

class [[maybe_unused]] PotentialCorrelation : public Analysis
{
  private:
    std::unique_ptr<pairpotential::NewCoulombGalore>
        coulomb; 
    Equidistant2DTable<double, Average<double>> mean_correlation; 
    std::string filename;            
    std::pair<double, double> range; 
    double dr = 0;                   
    unsigned int calculations_per_sample_event = 1;
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    PotentialCorrelation(const json& j, const Space& spc);
};

class ElectricPotential : public Analysis
{
  public:
    enum class Policies
    {
        FIXED,
        RANDOM_WALK,
        RANDOM_WALK_NO_OVERLAP,
        INVALID
    };

  private:
    unsigned int max_overlap_trials = 100; 
    double histogram_resolution = 0.01;    
    unsigned int calculations_per_sample_event = 1;
    std::string file_prefix; 

    struct Target
    {
        Point position;                 
        Average<double> mean_potential; 
        std::unique_ptr<SparseHistogram<double>>
            potential_histogram; 
    };

    std::vector<Target> targets; 
    Policies policy;             
    std::unique_ptr<pairpotential::NewCoulombGalore>
        coulomb; 
    Average<double>
        mean_potential_correlation; 
    SparseHistogram<double>
        potential_correlation_histogram; 
    void getTargets(const json& j);      
    void setPolicy(const json& j);       
    double
    calcPotentialOnTarget(const Target& target); 
    [[nodiscard]] bool overlapWithParticles(
        const Point& position) const;  
    std::function<void()> applyPolicy; 
    json output_information;           
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    ElectricPotential(const json& j, const Space& spc);
};

NLOHMANN_JSON_SERIALIZE_ENUM(ElectricPotential::Policies,
                             {
                                 {ElectricPotential::Policies::INVALID, nullptr},
                                 {ElectricPotential::Policies::FIXED, "fixed"},
                                 {ElectricPotential::Policies::RANDOM_WALK_NO_OVERLAP,
                                  "no_overlap"},
                                 {ElectricPotential::Policies::RANDOM_WALK, "random_walk"},
                             })

class AtomProfile : public Analysis
{
    Equidistant2DTable<double, double> table;         
    double dr;                                        
    std::vector<std::string> names;                   
    std::set<AtomData::index_type> atom_id_selection; 
    std::string file;                                 
    Point origin = {0.0, 0.0, 0.0};
    Eigen::Vector3i dir = {1, 1, 1};
    bool count_charge = false;
    int center_of_mass_atom_id = -1; // center at COM of id_com atoms?

    [[nodiscard]] double distanceToOrigin(const Point& position) const;
    void _from_json(const json& j) override;
    void _to_json(json& j) const override;
    void _to_disk() override;
    void _sample() override;

  public:
    AtomProfile(const json& j, const Space& spc);
};

class SlicedDensity : public Analysis
{
    double dz = 0.0;                         
    Table2D<double, unsigned int> histogram; // N(z)
    std::vector<std::string> atom_names;
    std::vector<AtomData::index_type> atom_ids; 
    std::string file;
    int center_of_mass_atom_id = -1; // center at COM of id_com atoms?

    void _from_json(const json& j) override;
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    SlicedDensity(const json& j, const Space& spc);
};

class Density : public Analysis
{
  protected:
    using Table = Equidistant2DTable<unsigned int, double>;
    using id_type = size_t;
    std::map<id_type, Average<double>> mean_density;
    std::map<id_type, std::string_view> names; // id <-> name database
    void _to_disk() override;
    void _sample() override;
    void _to_json(json& j) const override;
    static void writeTable(std::string_view name, Table& table);

  private:
    [[nodiscard]] virtual std::map<id_type, int> count() const = 0;
    std::map<MoleculeData::index_type, Table> probability_density;
    Average<double> mean_cubic_root_of_volume;
    Average<double> mean_volume;
    Average<double> mean_inverse_volume;
    double updateVolumeStatistics();

  public:
    template <RequireNamedElements Range>
    Density(const Space& spc, const Range& atoms_or_molecules, const std::string_view name)
        : Analysis(spc, name)
    {
        for (const auto& data : atoms_or_molecules) {
            names[data.id()] = data.name;
            probability_density[data.id()].setResolution(1, 0);
        }
    }
};

class MoleculeDensity final : public Density
{
  private:
    [[nodiscard]] std::map<id_type, int> count() const override;

  public:
    MoleculeDensity(const json& j, const Space& spc);
};

class AtomDensity final : public Density
{
  private:
    std::map<id_type, Table> atomswap_probability_density;
    void _sample() override;
    void _to_disk() override;
    [[nodiscard]] std::map<id_type, int> count() const override;

  public:
    AtomDensity(const json& j, const Space& spc);
};

std::function<double(const Group&, const Group&)>
createGroupGroupProperty(const json& j, const Space& spc, Energy::Hamiltonian& hamiltonian);

std::function<bool(double)> createValueFilter(const std::string& name, const json& j,
                                              bool throw_on_error);

class PairMatrixAnalysis : public Analysis
{
  private:
    std::unique_ptr<std::ostream> matrix_stream; 
    std::string filename;                        
    virtual void setPairMatrix() = 0;            
    void _to_json(json& j) const override;
    void _from_json(const json& j) override;
    void _sample() override;
    void _to_disk() override;

  protected:
    const Space& spc;
    std::function<bool(double)> value_filter; 
    Eigen::SparseMatrix<double> pair_matrix;  

  public:
    PairMatrixAnalysis(const json& j, const Space& spc);
};

class GroupMatrixAnalysis : public PairMatrixAnalysis
{
  private:
    std::function<double(const Group&, const Group&)>
        property;                   
    std::vector<int> group_indices; 
    void setPairMatrix() override;  
  public:
    GroupMatrixAnalysis(const json& j, const Space& spc, Energy::Hamiltonian& pot);
};

class ChargeFluctuations : public Analysis
{
  private:
    typename decltype(Faunus::molecules)::const_iterator mol_iter; 
    using AtomHistogram = std::map<AtomData::index_type, int>; 
    std::vector<AtomHistogram> atom_histograms;          
    std::vector<AverageStdev<double>> atom_mean_charges; 
    std::string filename;                                
    bool verbose = true;                                 

    ParticleVector averageChargeParticles(const Space::GroupType& group);
    [[nodiscard]] std::vector<double> getMeanCharges() const;
    [[nodiscard]] std::vector<double> getChargeStandardDeviation() const;
    [[nodiscard]] std::vector<std::string> getPredominantParticleNames() const;

    void _sample() override;
    void _to_json(json& j) const override;
    void _to_disk() override;

  public:
    ChargeFluctuations(const json& j, const Space& spc);
}; // Fluctuations of atomic charges

class Multipole : public Analysis
{
    struct Data
    {
        Average<double> charge;
        Average<double> charge_squared;
        Average<double> dipole_moment;
        Average<double> dipole_moment_squared;
    }; // Average sample moment for a molecule

    std::map<MoleculeData::index_type, Data>
        average_moments; 
    void _sample() override;
    void _to_json(json& j) const override;

  public:
    Multipole(const json& j, const Space& spc);
};

class SystemEnergy : public Analysis
{
  private:
    Energy::Hamiltonian& hamiltonian;
    std::string file_name;
    std::string separator;
    std::unique_ptr<std::ostream> output_stream;
    [[nodiscard]] std::vector<double> calculateEnergies() const;
    Average<double> mean_energy;
    Average<double> mean_squared_energy;
    Table2D<double, double> energy_histogram; // Density histograms
    double initial_energy = 0.0;
    double minimum_energy = std::numeric_limits<double>::infinity(); 
    bool dump_minimum_energy_configuration = false; 

    bool updateMinimumEnergy(double current_energy);
    void createOutputStream();
    [[maybe_unused]] void normalize();
    void _sample() override;
    void _to_json(json& j) const override;
    void _from_json(const json& j) override;
    void _to_disk() override;

  public:
    SystemEnergy(const json& j, const Space& spc, Energy::Hamiltonian& hamiltonian);
};

class SanityCheck : public Analysis
{
  private:
    const double mass_center_tolerance = 1.0e-3;
    void _sample() override;
    void checkGroupsCoverParticles(); 
    void
    checkMassCenter(const Space::GroupType& group); 
    void checkWithinContainer(
        const Space::GroupType& group); 

  public:
    SanityCheck(const json& j, const Space& spc);
};

class SaveState : public Analysis
{
  private:
    std::function<void(const std::string&)> writeFunc = nullptr;
    bool save_random_number_generator_state = false;
    bool use_numbered_files = true;
    bool convert_hexagonal_prism_to_cuboid = false;
    std::string filename;

    void _to_json(json& j) const override;
    void _sample() override;
    static void saveAsCuboid(const std::string& filename, const Space& spc,
                             StructureFileWriter& writer);
    void saveJsonStateFile(const std::string& filename, const Space& spc) const;
    void saveBinaryJsonStateFile(const std::string& filename, const Space& spc) const;

  public:
    SaveState(json j, const Space& spc);
    ~SaveState() override;
    void setWriteFunction(const Space& spc);
};

class PairFunction : public Analysis
{
  protected:
    using index_type = size_t;
    int dimensions = 3;                           
    index_type id1 = 0;                           
    index_type id2 = 0;                           
    double dr = 0;                                
    std::string name1;                            
    std::string name2;                            
    std::string file;                             
    Equidistant2DTable<double, double> histogram; 
    Average<double> mean_volume;                  
    Eigen::Vector3i slicedir = {0, 0, 0};
    double thickness = 0;

  private:
    void _from_json(const json& j) override;
    void _to_json(json& j) const override;
    void _to_disk() override;
    [[nodiscard]] double volumeElement(double r) const;

  public:
    PairFunction(const Space& spc, const json& j, std::string_view name);
};

class AtomRDF : public PairFunction
{
  private:
    void _sample() override;
    void sampleDistance(const Particle& particle1, const Particle& particle2);
    void sampleDifferent(); 
    void sampleIdentical(); 

  public:
    AtomRDF(const json& j, const Space& spc);
};

class MoleculeRDF : public PairFunction
{
  private:
    void _sample() override;
    void sampleDistance(const Group& group_i, const Group& group_j);
    void sampleDifferent(); 
    void sampleIdentical(); 

  public:
    MoleculeRDF(const json& j, const Space& spc);
};

class PairAngleFunction : public PairFunction
{
  private:
    std::string correlation_filename;

  protected:
    Equidistant2DTable<double, Average<double>> average_correlation_vs_distance;

  private:
    void _from_json(const json& j) override;
    void _to_disk() override;

  public:
    PairAngleFunction(const Space& spc, const json& j, const std::string& name);
};

class AtomDipDipCorr : public PairAngleFunction
{
    void _sample() override;

  public:
    AtomDipDipCorr(const json& j, const Space& spc);
};

class XTCtraj : public Analysis
{
  private:
    using index_type = MoleculeData::index_type;
    std::vector<index_type> group_ids;      
    std::vector<std::size_t> group_indices; 
    std::unique_ptr<XTCWriter> writer;
    void _to_json(json& j) const override;
    void _sample() override;
    XTCtraj(const Space& spc, const std::string& filename,
            const std::vector<std::string>& molecule_names);

  public:
    XTCtraj(const json& j, const Space& spc);
};

class VirtualVolumeMove : public PerturbationAnalysis
{
    Geometry::VolumeMethod volume_scaling_method = Geometry::VolumeMethod::ISOTROPIC;
    double volume_displacement = 0.0;
    void _sample() override;
    void _from_json(const json& j) override;
    void _to_json(json& j) const override;
    void sanityCheck(double old_energy);
    void writeToFileStream(const Point& scale, double energy_change) const;

  public:
    VirtualVolumeMove(const json& j, Space& spc, Energy::EnergyTerm& pot);
};

class MolecularConformationID : public Analysis
{
    MoleculeData::index_type molid;        
    std::map<int, unsigned int> histogram; 
    void _sample() override;
    void _to_json(json& j) const override;

  public:
    MolecularConformationID(const json& j, const Space& spc);
};

class VirtualTranslate : public PerturbationAnalysis
{
    MoleculeData::index_type molid; 
    Point perturbation_direction = {0.0, 0.0, 1.0};
    double perturbation_distance = 0.0;

    void _sample() override;
    void _from_json(const json& j) override;
    void _to_json(json& j) const override;
    double momentarilyPerturb(Space::GroupType& group);
    void writeToFileStream(double energy_change) const;

  public:
    VirtualTranslate(const json& j, Space& spc, Energy::EnergyTerm& pot);
};

class MultipoleDistribution : public Analysis
{
    struct Data
    {
        using average_type = Average<double>;
        average_type exact; 
        average_type ion_ion;
        average_type ion_dipole;
        average_type ion_quadrupole;
        average_type dipole_dipole;
        average_type dipole_dipole_correlation;
    };

    std::map<MoleculeData::index_type, Data> mean_energy; 
    std::vector<std::string> names;                       
    std::vector<MoleculeData::index_type> ids;            
    std::string filename;                                 
    double dr = 0.0;                                      

    void _sample() override;
    void _to_json(json& j) const override;
    void _to_disk() override;
    void sampleGroupGroup(const Space::GroupType& group1, const Space::GroupType& group2);
    [[nodiscard]] double groupGroupExactEnergy(
        const Space::GroupType& group1,
        const Space::GroupType& group2) const; //<! exact ion-ion energy between particles

  public:
    MultipoleDistribution(const json& j, const Space& spc);
}; // end of multipole distribution

class ScatteringFunction : public Analysis
{
  private:
    enum class Schemes
    {
        DEBYE,
        EXPLICIT_PBC,
        EXPLICIT_IPBC
    }; // three different schemes
    Schemes scheme = Schemes::DEBYE;
    bool mass_center_scattering;          
    bool save_after_sample = false;       
    std::string filename;                 
    std::vector<Scatter::Scatterer> scatter_positions; 
    std::vector<MoleculeData::index_type> molecule_ids; 
    std::vector<std::string> molecule_names; 
    using Tformfactor = Scatter::FormFactorAtomicConstant<double>;

    std::unique_ptr<Scatter::DebyeFormula<Tformfactor>> debye;
    std::unique_ptr<Scatter::StructureFactorPBC<Tformfactor>> explicit_average_pbc;
    std::unique_ptr<Scatter::StructureFactorIPBC<Tformfactor, float>> explicit_average_ipbc;
    void _sample() override;
    void _to_disk() override;
    void _to_json(json& j) const override;

  public:
    ScatteringFunction(const json& j, const Space& spc);
};

/*
 * @brief Sample and save gyration eigenvalues of all particles having the same id
 */
class AtomInertia : public Analysis
{
  private:
    std::string filename;
    AtomData::index_type atom_id;
    std::ofstream output_stream;

    Point compute();
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    AtomInertia(const json& j, const Space& spc);
};

class InertiaTensor : public Analysis
{
  private:
    std::string filename;                 
    std::unique_ptr<std::ostream> stream; 
    MoleculeData::index_type group_index; 
    std::vector<size_t> particle_range;   
    [[nodiscard]] std::pair<Point, Point>
    compute() const; 
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    InertiaTensor(const json& j, const Space& spc);
};

/*
 * @brief Sample and save charge, dipole, and quadrupole moments for a range of indexes within a
 * molecule
 */
class MultipoleMoments : public Analysis
{
  private:
    std::string filename;
    std::vector<AtomData::index_type> particle_range; 
    MoleculeData::index_type group_index;
    bool use_molecular_mass_center =
        true; 
    std::ofstream output_stream;

    struct Data
    {
        double charge = 0.0;                // total charge
        Point dipole_moment{0.0, 0.0, 0.0}; // dipole vector
        Point eivals, eivec, center;        // quadrupole eigenvalues and major axis
    } __attribute__((aligned(128)));

    Data calculateMultipoleMoment() const;
    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    MultipoleMoments(const json& j, const Space& spc);
};

class PolymerShape : public Analysis
{
    struct AverageData
    {
        using average_type = Average<double>;
        average_type gyration_radius_squared;
        average_type gyration_radius;
        average_type end_to_end_squared;
        average_type shape_factor_squared;
        average_type aspherity;
        average_type acylindricity;
        average_type relative_shape_anisotropy;
    }; 

    AverageData data;               
    MoleculeData::index_type molid; 
    std::unique_ptr<std::ostream> tensor_output_stream = nullptr;
    Equidistant2DTable<double, unsigned int> gyration_radius_histogram;

    void _to_json(json& j) const override;
    void _sample() override;
    void _to_disk() override;

  public:
    PolymerShape(const json& j, const Space& spc);
};

class QRtraj : public Analysis
{
  protected:
    std::function<void()> write_to_file;            
    std::unique_ptr<std::ostream> stream = nullptr; 

  private:
    std::string filename;    
    void _sample() override; 
    void _to_json(json& j) const override;
    void _to_disk() override;

  public:
    QRtraj(const json& j, const Space& spc, const std::string& name = "qrtraj");
};

class PatchySpheroCylinderTrajectory : public QRtraj
{
  public:
    PatchySpheroCylinderTrajectory(const json& j, const Space& spc);
};

class AreaSamplingPolicy;

class SASAAnalysis : public Analysis
{
    using index_type = Faunus::AtomData::index_type;
    using count_type = size_t;
    using table_type = Equidistant2DTable<double, count_type>;

  public:
    enum class Policies
    {
        ATOMIC,
        MOLECULAR,
        ATOMS_IN_MOLECULE,
        INVALID
    };

  private:
    struct Averages
    {
        using average_type = Average<double>;
        average_type area;
        average_type area_squared;
    }; 

    Averages average_data; 

    std::unique_ptr<std::ostream> output_stream; 

    double probe_radius; 
    int slices_per_atom; 

    std::string filename;                         
    table_type sasa_histogram;                    
    std::unique_ptr<Faunus::SASA::SASABase> sasa; 
    std::unique_ptr<AreaSamplingPolicy> policy; 

    void _to_json(json& j) const override;
    void _from_json(const json& input) override;
    void _to_disk() override;
    void _sample() override;

    void setPolicy(Policies);
    void takeSample(double area);

    friend class AreaSamplingPolicy;

  public:
    SASAAnalysis(const json& j, const Space& spc);
    SASAAnalysis(double probe_radius, int slices_per_atom, double resolution, Policies policy,
                 const Space& spc);
};

NLOHMANN_JSON_SERIALIZE_ENUM(SASAAnalysis::Policies,
                             {{SASAAnalysis::Policies::ATOMIC, "atomic"},
                              {SASAAnalysis::Policies::MOLECULAR, "molecular"},
                              {SASAAnalysis::Policies::ATOMS_IN_MOLECULE, "atoms_in_molecule"},
                              {SASAAnalysis::Policies::INVALID, nullptr}})

class AreaSamplingPolicy
{
  protected:
    template <typename TBegin, typename TEnd>
    void sampleIndividualSASA(TBegin first, TEnd last, SASAAnalysis& analysis);
    template <typename TBegin, typename TEnd>
    void sampleTotalSASA(TBegin first, TEnd last, SASAAnalysis& analysis);

  public:
    virtual void sample(const Space& spc, SASAAnalysis& analysis) = 0;
    virtual void to_json(json& input) const = 0;
    virtual void from_json(const json& input) = 0;

    inline AreaSamplingPolicy() = default;
    inline virtual ~AreaSamplingPolicy() = default;
};

class AtomicPolicy final : public AreaSamplingPolicy
{

    AtomData::index_type atom_id; 
    std::string atom_name;        

    void sample(const Space& spc, SASAAnalysis& analysis) override;
    void to_json(json& input) const override;
    void from_json(const json& input) override;

  public:
    inline AtomicPolicy() = default;
};

class MolecularPolicy final : public AreaSamplingPolicy
{

    MoleculeData::index_type molecule_id; 
    std::string molecule_name;            

    void sample(const Space& spc, SASAAnalysis& analysis) override;
    void to_json(json& input) const override;
    void from_json(const json& input) override;

  public:
    inline MolecularPolicy() = default;
};

class AtomsInMoleculePolicy final : public AreaSamplingPolicy
{

    MoleculeData::index_type molecule_id; 
    std::set<size_t> selected_indices;    
    std::string molecule_name;            
    std::set<std::string> atom_names;     

    void sample(const Space& spc, SASAAnalysis& analysis) override;
    void to_json(json& input) const override;
    void from_json(const json& input) override;

  public:
    inline AtomsInMoleculePolicy() = default;
};

template <class T, class Enable = void> struct [[maybe_unused]] _analyse
{
    void sample(T&) { std::cout << "not a dipole!" << std::endl; } 
}; // primary template

template <class T>
struct [[maybe_unused]] _analyse<T,
                                 typename std::enable_if<std::is_base_of<Dipole, T>::value>::type>
{
    void sample(T&) { std::cout << "dipole!" << std::endl; } 
}; // specialized template

class SavePenaltyEnergy : public Analysis
{
  private:
    const std::string filename; 
    int filenumber = 0;         
    std::shared_ptr<Energy::Penalty> penalty_energy;
    void _sample() override;

  public:
    SavePenaltyEnergy(const json& j, const Space& spc, const Energy::Hamiltonian& pot);
};

class Voronota : public Analysis
{
  public:
    enum class Modes
    {
        FULL,
        INTERCHAIN,
        UPDATEABLE,
        INVALID
    };

  private:
    struct Averages
    {
        Average<double> area_of_contacts;
        Average<double> volume_inside_sas;
        Average<double> area_of_sas;
        Average<double> area_of_sas_per_chain;
    }; 

    Modes mode = Modes::FULL;                    
    double probe_radius;                         
    bool use_pbc = false;                        
    Averages average_data;                       
    std::unique_ptr<std::ostream> output_stream; 
    std::string filename;                        

    class Impl;
    std::unique_ptr<Impl> pimpl;

    void _to_json(json& json_output) const override;
    void _from_json(const json& input) override;
    void _to_disk() override;
    void _sample() override;

  public:
    Voronota(const json& j, const Space& spc);
    Voronota(Modes mode, double probe_radius, const Space& spc);
    ~Voronota();
};

} // namespace Faunus::analysis