#if defined(DEBUG) && not defined(RELEASE_BRANCH)
#include <fstream>
#include <filesystem>
#endif

#include <calculate_integral.hpp>
#include <algorithm/operate_on_merged_set_intervals.hpp>

inline uint32_t scan_interval_boarder(double l, double r, stl_vector_mp<double>& parallel_line)
{
    std::sort(parallel_line.begin(), parallel_line.end());
    auto unique_end_iter = std::unique(parallel_line.begin(), parallel_line.end(), sorted_double_equal);
    parallel_line.erase(unique_end_iter, parallel_line.end());

    auto begin_iter =
        std::find_if(parallel_line.begin(), parallel_line.end(), [&](double x) { return x - l > strict_epsilon; });
    auto end_iter = std::find_if(parallel_line.begin(), parallel_line.end(), [&](double x) { return r - x <= strict_epsilon; });
    return std::distance(begin_iter, end_iter);
}

struct chain_vertex_property_t {
    uint64_t total_count : 32;
    uint64_t neg_delta_u : 8;
    uint64_t pos_delta_u : 8;
    uint64_t neg_delta_v : 8;
    uint64_t pos_delta_v : 8;

    inline chain_vertex_property_t& operator+=(const chain_vertex_property_t& x) noexcept
    {
        this->total_count += x.total_count;
        this->neg_delta_u += x.neg_delta_u;
        this->pos_delta_u += x.pos_delta_u;
        this->neg_delta_v += x.neg_delta_v;
        this->pos_delta_v += x.pos_delta_v;
        return *this;
    }

    inline chain_vertex_property_t operator~() noexcept
    {
        chain_vertex_property_t res{};
        res.total_count = this->total_count;
        res.neg_delta_u = this->pos_delta_u;
        res.pos_delta_u = this->neg_delta_u;
        res.neg_delta_v = this->pos_delta_v;
        res.pos_delta_v = this->neg_delta_v;
        return res;
    }
};

void calculate_integral(uint8_t                                        q,
                        const baked_blobtree_t&                        tree,
                        const pcurve_relation_graph_t&                 pcurve_relation_graph,
                        stl_vector_mp<stl_vector_mp<Eigen::Vector4d>>& integral_points,
                        stl_vector_mp<stl_vector_mp<double>>&          surface_integral_weights,
                        stl_vector_mp<stl_vector_mp<double>>&          volume_integral_weights)
{
    using algorithm::operate_on_merged_set_intervals;

    const auto& subfaces      = tree.subfaces;
    const auto& subface_types = tree.subface_types;

    flat_hash_map_mp<Eigen::Vector2d, chain_vertex_property_t, chain_end_hasher> chain_ends{};
    stl_vector_mp<std::pair<Eigen::Vector2d, chain_vertex_property_t>>           subchain_ends{};
    stl_vector_mp<double>                                                        line_parallel_to_u{}, line_parallel_to_v{};
    stl_vector_mp<double>                                                        intersection_u{}, intersection_v{};
    stl_vector_mp<Eigen::Vector3d>                                               param_integral_points{};

    auto tensor_product_result = calculate_tensor_product_points(q);

    const auto& subface_nodes = pcurve_relation_graph.nodes<0>();
    integral_points.resize(subface_nodes.size());
    surface_integral_weights.resize(subface_nodes.size());
    volume_integral_weights.resize(subface_nodes.size());
    for (auto subface_node : subface_nodes) {
        // early exit: iff no active patch
        if (subface_node.edges<edge_bidirection::to_next>().empty()) continue;

        chain_ends.clear();
        line_parallel_to_u.clear();
        line_parallel_to_v.clear();
        intersection_u.clear();
        intersection_v.clear();
        param_integral_points.clear();

        const auto subface_index = subface_node.index();
        const auto subface_ptr   = subfaces[subface_index];
        const auto subface_type  = subface_types[subface_index];
        auto       map_param_to_point_with_weight =
            std::bind(internal::get_map_param_to_point_with_weight_ptr(subface_type), subface_ptr, std::placeholders::_1);
        const auto subface_range_u = internal::get_range_u(subface_type), subface_range_v = internal::get_range_v(subface_type);
        const auto u_range = subface_range_u.first, v_range = subface_range_v.first;

        auto& int_p                     = integral_points[subface_index];
        auto& sur_int_w                 = surface_integral_weights[subface_index];
        auto& vol_int_w                 = volume_integral_weights[subface_index];
        auto  transform_integral_points = [&] {
            std::for_each(param_integral_points.begin(), param_integral_points.end(), [&](const Eigen::Vector3d& p) {
                auto [world_p, surface_weight, volume_weight] = map_param_to_point_with_weight(Eigen::Vector2d{p.x(), p.y()});
                int_p.emplace_back(std::move(world_p));
                sur_int_w.emplace_back(p.z() * surface_weight);
                vol_int_w.emplace_back(p.z() * volume_weight);
            });
        };

        auto edges_to_chains  = subface_node.edges<edge_bidirection::to_prev>();
        auto edges_to_patches = subface_node.edges<edge_bidirection::to_next>();
        // if no chain on subface, then just return tensor-product points as result
        // only one condition is permitted: subface is self sealing, as sphere
        if (bool empty_chains = edges_to_chains.empty(), limited_range = u_range < infinity && v_range < infinity;
            empty_chains && limited_range) {
            // if this subface is permitted to be plane, then just take use of tensor product points
            int_p.reserve(q * q);
            sur_int_w.reserve(q * q);
            vol_int_w.reserve(q * q);
            std::for_each(tensor_product_result.begin(), tensor_product_result.end(), [&](const Eigen::Vector3d& p) {
                Eigen::Vector2d transform_p{p.x() * u_range, p.y() * v_range};
                auto [world_p, surface_weight, volume_weight] = map_param_to_point_with_weight(transform_p);
                int_p.emplace_back(std::move(world_p));
                sur_int_w.emplace_back(p.z() * surface_weight);
                vol_int_w.emplace_back(p.z() * volume_weight);
            });
            continue;
        } else if (empty_chains)
            assert(limited_range);

        Eigen::AlignedBox2d subface_aabb{};
        for (auto edge_to_chain : edges_to_chains)
            for (const auto& [_, aabb] : edge_to_chain.property().subchains) subface_aabb.extend(aabb);
        // EDIT: change aabb bound to uv range if they exist, so we can include any valid interval
        if (u_range < infinity) {
            subface_aabb.min().x() = .0;
            subface_aabb.max().x() = u_range;
        }
        if (v_range < infinity) {
            subface_aabb.min().y() = .0;
            subface_aabb.max().y() = v_range;
        }

        // on boarder u = x
        auto scan_u_boarder_interval = [&](double x) {
            return u_line_subface_intersection(pcurve_relation_graph,
                                               subface_index,
                                               x,
                                               subface_aabb.min().y(),
                                               subface_aabb.max().y(),
                                               true,
                                               true);
        };
        // on boarder v = x
        auto scan_v_boarder_interval = [&](double x) {
            return v_line_subface_intersection(pcurve_relation_graph,
                                               subface_index,
                                               x,
                                               subface_aabb.min().x(),
                                               subface_aabb.max().x(),
                                               true,
                                               true);
        };

        auto expand_integral_on_u_ = [&](auto i, double l, double r, auto b0_iter, auto b1_iter, auto int_x, auto int_w) {
            auto x = gl_integrator.x_interval(q, i, l, r);
            if (x - b0_iter[1] > strict_epsilon) b0_iter += 2;
            if (x - b1_iter[1] > strict_epsilon) b1_iter += 2;
            expand_integral_on_u(q,
                                 pcurve_relation_graph,
                                 subface_index,
                                 x,
                                 subface_aabb.min().x(),
                                 subface_aabb.max().x(),
                                 x - b0_iter[0] > strict_epsilon,
                                 x - b1_iter[0] > strict_epsilon,
                                 gl_integrator.w_interval(q, i, l, r),
                                 int_x,
                                 int_w,
                                 param_integral_points);
        };
        auto expand_integral_on_v_ = [&](auto i, double l, double r, auto b0_iter, auto b1_iter, auto int_x, auto int_w) {
            auto x = gl_integrator.x_interval(q, i, l, r);
            if (x - b0_iter[1] > strict_epsilon) b0_iter += 2;
            if (x - b1_iter[1] > strict_epsilon) b1_iter += 2;
            expand_integral_on_v(q,
                                 pcurve_relation_graph,
                                 subface_index,
                                 x,
                                 subface_aabb.min().y(),
                                 subface_aabb.max().y(),
                                 x - b0_iter[0] > strict_epsilon,
                                 x - b1_iter[0] > strict_epsilon,
                                 gl_integrator.w_interval(q, i, l, r),
                                 int_x,
                                 int_w,
                                 param_integral_points);
        };
        auto expand_integral_on_u_by_gl = [&](auto i, double l, double r, auto b0_iter, auto b1_iter) {
            expand_integral_on_u_(i, l, r, b0_iter, b1_iter, gl_integrator.x_interval, gl_integrator.w_interval);
        };
        auto expand_integral_on_u_by_ts = [&](auto i, double l, double r, auto b0_iter, auto b1_iter) {
            expand_integral_on_u_(i, l, r, b0_iter, b1_iter, ts_integrator.x_interval, ts_integrator.w_interval);
        };
        auto expand_integral_on_v_by_gl = [&](auto i, double l, double r, auto b0_iter, auto b1_iter) {
            expand_integral_on_v_(i, l, r, b0_iter, b1_iter, gl_integrator.x_interval, gl_integrator.w_interval);
        };
        auto expand_integral_on_v_by_ts = [&](auto i, double l, double r, auto b0_iter, auto b1_iter) {
            expand_integral_on_v_(i, l, r, b0_iter, b1_iter, ts_integrator.x_interval, ts_integrator.w_interval);
        };

        for (auto edge_to_chain : edges_to_chains) {
            const auto  chain_index = edge_to_chain.to().index();
            const auto& subchains   = edge_to_chain.property().subchains;
            subchain_ends.clear();
            for (uint32_t subchain_index = 0; subchain_index < subchains.size(); subchain_index++) {
                const auto& vertices = subchains[subchain_index].vertices;

                auto vertex_iter      = vertices.begin();
                auto prev_vertex_iter = vertex_iter;
                vertex_iter++;
                chain_vertex_property_t vertex_property{};
                auto                    update_vertex_property = [&] {
                    auto prev_to_cur_delta      = *vertex_iter - *prev_vertex_iter;
                    vertex_property             = {};
                    vertex_property.total_count = 1;

                    // condition below: negative, near parallel, positive
                    if (prev_to_cur_delta.x() < -strict_epsilon)
                        vertex_property.neg_delta_u = 1;
                    else if (prev_to_cur_delta.x() >= strict_epsilon)
                        vertex_property.pos_delta_u = 1;
                    else
                        line_parallel_to_v.emplace_back(vertex_iter->x());

                    if (prev_to_cur_delta.y() < -strict_epsilon)
                        vertex_property.neg_delta_v = 1;
                    else if (prev_to_cur_delta.y() >= strict_epsilon)
                        vertex_property.pos_delta_v = 1;
                    else
                        line_parallel_to_u.emplace_back(vertex_iter->y());
                };
                update_vertex_property();
                subchain_ends.emplace_back(*prev_vertex_iter, vertex_property);
                prev_vertex_iter++, vertex_iter++;
                auto prev_vertex_property = vertex_property;
                for (; prev_vertex_iter != vertices.end() - 1; prev_vertex_iter++, vertex_iter++) {
                    update_vertex_property();
                    if ((prev_vertex_property.neg_delta_u && vertex_property.pos_delta_u) //
                        || (prev_vertex_property.pos_delta_u && vertex_property.neg_delta_u))
                        line_parallel_to_v.emplace_back(vertex_iter->x());
                    if ((prev_vertex_property.neg_delta_v && vertex_property.pos_delta_v) //
                        || (prev_vertex_property.pos_delta_v && vertex_property.neg_delta_v))
                        line_parallel_to_u.emplace_back(vertex_iter->y());

                    prev_vertex_property = vertex_property;
                }
                subchain_ends.emplace_back(*prev_vertex_iter, ~vertex_property);
            }

            // handle subchain end, and insert ends into chain ends
            auto subchain_end_iter                = subchain_ends.begin();
            chain_ends[subchain_end_iter->first] += subchain_end_iter->second;
            subchain_end_iter++;
            for (; subchain_end_iter != subchain_ends.end() - 1; subchain_end_iter += 2) {
                // HINT: the vertex of iter, iter + 1 should be the same except for periodical remapping
                const auto& p         = subchain_end_iter->first;
                auto&       property  = subchain_end_iter->second;
                property             += (subchain_end_iter + 1)->second;

                if (property.neg_delta_u >= 2 || property.pos_delta_u >= 2) line_parallel_to_v.emplace_back(p.x());
                if (property.neg_delta_v >= 2 || property.pos_delta_v >= 2) line_parallel_to_u.emplace_back(p.y());
            }
            // CAUTION: avoid ring
            if (subchain_ends.front().first != subchain_ends.back().first)
                chain_ends[subchain_end_iter->first] += subchain_end_iter->second;
        }

        // handle each chain end, try to figure out if there is intersection point or near parallel tagent
        for (const auto& [p, property] : chain_ends) {
            if (property.total_count > 1) {
                intersection_u.emplace_back(p.x());
                intersection_v.emplace_back(p.y());
            }

            if (property.neg_delta_u >= 2 || property.pos_delta_u >= 2) line_parallel_to_v.emplace_back(p.x());
            if (property.neg_delta_v >= 2 || property.pos_delta_v >= 2) line_parallel_to_u.emplace_back(p.y());
        }

        const auto line_parallel_to_u_size =
            scan_interval_boarder(subface_aabb.min().y(), subface_aabb.max().y(), line_parallel_to_u);
        const auto line_parallel_to_v_size =
            scan_interval_boarder(subface_aabb.min().x(), subface_aabb.max().x(), line_parallel_to_v);

#if defined(DEBUG) && not defined(RELEASE_BRANCH)
        auto path = std::filesystem::current_path().parent_path().parent_path().parent_path().parent_path()
                    / "surface_integral_v2" / "plot";
        // append to given subface's meta data
        std::ofstream file(path / std::string("subface_" + std::to_string(subface_index) + ".meta"),
                           std::ios::out | std::ios::app);

        file << std::scientific << std::setprecision(16) << subface_aabb.min().x() << " ";
        file << std::scientific << std::setprecision(16) << subface_aabb.min().y() << " ";
        file << std::scientific << std::setprecision(16) << subface_aabb.max().x() << " ";
        file << std::scientific << std::setprecision(16) << subface_aabb.max().y() << " ";
        file << "\n\n";

        // output parallel line u = x
        for (const auto u : line_parallel_to_v) file << std::scientific << std::setprecision(16) << u << " ";
        file << "\n";
        for (const auto u : intersection_u) file << std::scientific << std::setprecision(16) << u << " ";
        file << "\n";
        // output parallel line v = x
        for (const auto v : line_parallel_to_u) file << std::scientific << std::setprecision(16) << v << " ";
        file << "\n";
        for (const auto v : intersection_v) file << std::scientific << std::setprecision(16) << v << " ";
        file << "\n";

        file.close();
#endif

        // always expand on direction of less parallel line, since this will generate less interval for integral
        // HINT: if less line parallel to direction u, then the integral should expand on direction v
        // iff either side of the interval is near parallel line, then corresponding half of the interval should use tanh-sinh
        // integrator otherwise we will always use gauss-legendre integrator
        if (line_parallel_to_v_size <= line_parallel_to_u_size) {
            intersection_u.emplace_back(subface_aabb.min().x());
            intersection_u.emplace_back(subface_aabb.max().x());
            std::sort(intersection_u.begin(), intersection_u.end());
            auto unique_end              = std::unique(intersection_u.begin(), intersection_u.end(), sorted_double_equal);
            auto start_boarder_intervals = scan_v_boarder_interval(subface_aabb.min().y());
            auto end_boarder_intervals   = scan_v_boarder_interval(subface_aabb.max().y());
            auto b0_iter = start_boarder_intervals.begin(), b1_iter = end_boarder_intervals.begin();
            operate_on_merged_set_intervals(line_parallel_to_v.begin(),
                                            line_parallel_to_v.end(),
                                            intersection_u.begin(),
                                            unique_end,
                                            double_less,
                                            [&](double l, double r, bool l_exact_gl, bool r_exact_gl) {
                                                if (l_exact_gl)
                                                    for (auto i = 0; i < q / 2; ++i)
                                                        expand_integral_on_v_by_gl(i, l, r, b0_iter, b1_iter);
                                                else
                                                    for (auto i = 0; i < q / 2; ++i)
                                                        expand_integral_on_v_by_ts(i, l, r, b0_iter, b1_iter);

                                                if (r_exact_gl)
                                                    for (auto i = q / 2; i < q; ++i)
                                                        expand_integral_on_v_by_gl(i, l, r, b0_iter, b1_iter);
                                                else
                                                    for (auto i = q / 2; i < q; ++i)
                                                        expand_integral_on_v_by_ts(i, l, r, b0_iter, b1_iter);
                                            });
            // integral_plan.axis_to_expand = axis::v;
        } else {
            intersection_v.emplace_back(subface_aabb.min().y());
            intersection_v.emplace_back(subface_aabb.max().y());
            std::sort(intersection_v.begin(), intersection_v.end());
            auto unique_end              = std::unique(intersection_v.begin(), intersection_v.end(), sorted_double_equal);
            auto start_boarder_intervals = scan_u_boarder_interval(subface_aabb.min().x());
            auto end_boarder_intervals   = scan_u_boarder_interval(subface_aabb.max().x());
            auto b0_iter = start_boarder_intervals.begin(), b1_iter = end_boarder_intervals.begin();
            operate_on_merged_set_intervals(line_parallel_to_u.begin(),
                                            line_parallel_to_u.end(),
                                            intersection_v.begin(),
                                            unique_end,
                                            double_less,
                                            [&](double l, double r, bool l_exact_gl, bool r_exact_gl) {
                                                if (l_exact_gl)
                                                    for (auto i = 0; i < q / 2; ++i)
                                                        expand_integral_on_u_by_gl(i, l, r, b0_iter, b1_iter);
                                                else
                                                    for (auto i = 0; i < q / 2; ++i)
                                                        expand_integral_on_u_by_ts(i, l, r, b0_iter, b1_iter);

                                                if (r_exact_gl)
                                                    for (auto i = q / 2; i < q; ++i)
                                                        expand_integral_on_u_by_gl(i, l, r, b0_iter, b1_iter);
                                                else
                                                    for (auto i = q / 2; i < q; ++i)
                                                        expand_integral_on_u_by_ts(i, l, r, b0_iter, b1_iter);
                                            });
            // integral_plan.axis_to_expand = axis::u;
        }

#if defined(DEBUG) && not defined(RELEASE_BRANCH)
        // append to given subface's meta data
        file.open(path / std::string("subface_" + std::to_string(subface_index) + ".int2d"));

        for (const auto& integral_points : param_integral_points)
            file << integral_points[0] << " " << integral_points[1] << " " << integral_points[2] << "\n";

        file.close();
#endif

        // transform integral points in [u, v] to [x, y, z] with related weights
        transform_integral_points();

        int_p.shrink_to_fit();
        sur_int_w.shrink_to_fit();
        vol_int_w.shrink_to_fit();

#if defined(DEBUG) && not defined(RELEASE_BRANCH)
        // append to given subface's meta data
        file.open(path / std::string("subface_" + std::to_string(subface_index) + ".int3d"));

        for (const auto& integral_points : int_p)
            file << integral_points[0] << " " << integral_points[1] << " " << integral_points[2] << "\n";
        file << "\n";

        for (const auto& weight : sur_int_w) file << weight << "\n";
        file << "\n";

        for (const auto& weight : vol_int_w) file << weight << "\n";
        file << "\n";

        file.close();
#endif
    }
}