#pragma once

#include <base/subface.hpp>
#include <process.hpp>
#include <container/stl_alias.hpp>
#include <container/wrapper/flat_index_group.hpp>

namespace internal
{

/**
 * @brief Numerical integrator for parametric surfaces with trimming curves.
 *
 * This class computes integrals (e.g., area, mass, etc.) over trimmed parametric surfaces
 * by subdividing the parameter domain and applying Gaussian quadrature.
 *
 * The integrator does not own the input data; it holds const references to ensure zero-copy semantics.
 * Users must ensure that the lifetime of input data exceeds that of the integrator.
 */
class SI_API integrator_t
{
public:
    /**
     * @note This constructor does not copy the data; it stores const references.
     *       The caller is responsible for ensuring the validity of the referenced data
     *       throughout the lifetime of this integrator.
     */
    integrator_t(const stl_vector_mp<object_with_index_mapping<subface>>& surfaces,
                 const flat_hash_map_mp<uint32_t, parametric_plane_t>&    uv_planes);

    /// Default destructor
    ~integrator_t() = default;


    double calculate(int gauss_order,
                     double (*func)(double                 u,
                                    double                 v,
                                    const Eigen::Vector3d& p,
                                    const Eigen::Vector3d& dU,
                                    const Eigen::Vector3d& dV)) const;


    double calculate_one_subface(const subface&                            subface,
                                 const parametric_plane_t&                 param_plane,
                                 const stl_vector_mp<Eigen::AlignedBox2d>& chain_bboxes,
                                 int                                       gauss_order,
                                 double (*func)(double                 u,
                                                double                 v,
                                                const Eigen::Vector3d& p,
                                                const Eigen::Vector3d& dU,
                                                const Eigen::Vector3d& dV)) const;

private:
    /// Non-owning reference to the list of subfaces
    const stl_vector_mp<object_with_index_mapping<subface>>& m_subfaces;

    /// Non-owning reference to the map of parametric planes (ID -> parametric_plane_t)
    const flat_hash_map_mp<uint32_t, parametric_plane_t>& m_uv_planes;

    // Precomputed bounding boxes for each chain in each parametric plane
    flat_hash_map_mp<uint32_t, stl_vector_mp<Eigen::AlignedBox2d>> m_chain_bboxes_hash{};

    stl_vector_mp<double> compute_u_breaks(const parametric_plane_t& param_plane, double u_min, double u_max) const;

    void find_v_intersections_at_u(const subface&                            subface,
                                   const parametric_plane_t&                 param_plane,
                                   const stl_vector_mp<Eigen::AlignedBox2d>& chain_bboxes,
                                   double                                    u_val,
                                   double                                    v_min,
                                   double                                    v_max,
                                   stl_vector_mp<double>&                    intersections,
                                   stl_vector_mp<uint16_t>&                  intersected_chains) const;


    bool is_edge_inside_domain(const stl_vector_mp<stl_vector_mp<uint16_t>>& chain_group_indices,
                               uint16_t                                      chain_idx1,
                               uint16_t                                      chain_idx2) const;

    bool is_u_near_singularity(double u, double tol = 1e-6) const;


    void sort_and_unique_with_tol(stl_vector_mp<double>& vec, double epsilon = 1e-8) const;
};

double newton_method(const std::function<implicit_equation_intermediate(double)>& F,
                     double                                                       v_initial,
                     double                                                       tolerance      = 1e-8,
                     int                                                          max_iterations = 100);

double area_integrand(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)
{
    return 1.0;
}

// Integrate (1/3) * r · n dA over all subfaces
double volume_integrand(double u, double v, const Eigen::Vector3d& p, const Eigen::Vector3d& dU, const Eigen::Vector3d& dV)
{
    Eigen::Vector3d cross = dU.cross(dV);
    return (p.dot(cross)) / 3.0; // (1/3) * (x dydz + y dzdx + z dxdy)
};


} // namespace internal