use crate::HowzatError as Error;
use crate::dd::Umpire;
use crate::dd::index::RayIncidenceIndex as _;
use crate::dd::mode::{HalfspaceMode, UmpireHalfspaceMode};
use crate::dd::ray::{
    AdjacencyEdge, EdgeTarget, RayData as _, RayId, RayListHeads, RayPartition, RayPartitionOwned,
};
use crate::dd::state::{ConeDd, ConeEngine, ConeOutput};
use crate::dd::zero::{ZeroRepr, ZeroSet as _};
use calculo::num::{Num, Sign};
use hullabaloo::types::{ComputationStatus, Representation, RepresentationKind, Row, RowSet};

#[derive(Clone, Copy)]
enum DegeneracyCandidates<'a> {
    Slice(&'a [RayId]),
    BitSet(&'a RowSet),
}

impl<N: Num, R: Representation, U, ZR> ConeDd<N, R, U, ZR>
where
    U: Umpire<N, ZR>,
    ZR: crate::dd::mode::PreorderedBackend,
{
    pub fn run_to_completion(mut self) -> Result<ConeOutput<N, R, U, ZR>, Error> {
        loop {
            if self.state.core.recompute_row_order {
                let (ctx, umpire, strict) = (
                    &mut self.state.core.ctx,
                    &mut self.state.umpire,
                    &self.state.core._strict_inequality_set,
                );
                umpire.recompute_row_order_vector(ctx, strict);
                let priority = self.state.core.weakly_added_halfspaces.clone();
                self.state
                    .umpire
                    .bump_priority_rows(&mut self.state.core.ctx, &priority);
                <UmpireHalfspaceMode<N, U, ZR> as HalfspaceMode>::on_row_order_recomputed(
                    &mut self.state,
                );
                #[cfg(debug_assertions)]
                {
                    self.state.debug_assert_iteration_prefix();
                }
                self.state.core.recompute_row_order = false;
                self.state.core.edges_dirty = true;
            }
            if self.state.core.comp_status == ComputationStatus::RegionEmpty
                || self.state.core.comp_status == ComputationStatus::AllFound
            {
                break;
            }
            let active = self.state.core.ray_graph.active_len();
            let Some(hh) = self.state.umpire.choose_next_halfspace(
                &self.state.core.ctx,
                &self.state.core.weakly_added_halfspaces,
                self.state.core.iter_state.iteration,
                active,
            ) else {
                self.state.core.comp_status = ComputationStatus::AllFound;
                break;
            };
            <UmpireHalfspaceMode<N, U, ZR> as HalfspaceMode>::add_halfspace(
                &mut self.state,
                hh,
            )?;
        }
        Ok(ConeOutput { state: self.state })
    }
}

impl<N: Num, R: Representation, U, ZR> ConeEngine<N, R, U, ZR>
where
    U: Umpire<N, ZR>,
    ZR: ZeroRepr,
{
    pub(crate) fn update_ray_orders(&mut self) {
        let relaxed = self.core.options.relaxed_enumeration();
        let mut ids = std::mem::take(&mut self.core.active_id_scratch);
        self.core.ray_graph.copy_active_ids(&mut ids);

        let mut changed = false;
        for id in ids.iter().copied() {
            let Some(ray_data) = self.core.ray_graph.ray_data_mut(id) else {
                continue;
            };
            let old = ray_data.first_infeasible_row();
            self.umpire
                .update_first_infeasible_row(&self.core.ctx, ray_data, relaxed);
            changed |= old != ray_data.first_infeasible_row();
        }
        self.core.active_id_scratch = ids;
        if changed {
            self.core.edges_dirty = true;
        }
    }

    pub fn ray_count(&self) -> usize {
        self.core.ray_graph.active_len()
    }

    #[inline(always)]
    fn take_zero_set(&mut self) -> <ZR as ZeroRepr>::Set {
        self.core
            .zero_set_pool
            .pop()
            .unwrap_or_else(|| ZR::empty_set(self.row_count()))
    }

    #[inline(always)]
    fn recycle_zero_set(&mut self, zero_set: <ZR as ZeroRepr>::Set) {
        self.core.zero_set_pool.push(zero_set);
    }

    #[inline(always)]
    fn insert_added_halfspace_zero(&mut self, row: Row) {
        if let Some(id) = ZR::id_for_row(&self.core.ctx, row) {
            if !self.core.added_halfspaces_zero.contains(id) {
                self.core.added_halfspaces_zero.insert(id);
                self.core.added_halfspaces_zero_size += 1;
            }
        }
    }

    fn rebuild_added_halfspaces_zero(&mut self) {
        self.core.added_halfspaces_zero.clear();
        self.core.added_halfspaces_zero_size = 0;
        let equality_rows: Vec<Row> = self
            .core
            .equality_set
            .iter()
            .map(|row| row.as_index())
            .collect();
        for row in equality_rows {
            self.insert_added_halfspace_zero(row);
        }
        let rows: Vec<Row> = self
            .core
            .added_halfspaces
            .iter()
            .map(|row| row.as_index())
            .collect();
        for row in rows {
            self.insert_added_halfspace_zero(row);
        }
    }

    #[inline(always)]
    fn recycle_zero_set_from_ray_data(&mut self, ray_data: &mut U::RayData) {
        let mut set = ZR::empty_set(0);
        std::mem::swap(ray_data.zero_set_mut(), &mut set);
        self.recycle_zero_set(set);
    }

    fn position_of_row(&self, row: Row) -> Option<Row> {
        self.core
            .ctx
            .row_to_pos
            .get(row)
            .copied()
            .filter(|pos| *pos < self.row_count())
    }

    pub(crate) fn first_infeasible_position(&self, ray: &U::RayData) -> Option<Row> {
        ray.first_infeasible_row()
            .and_then(|row| self.position_of_row(row))
    }

    fn first_infeasible_position_or_m(&self, ray: &U::RayData) -> Row {
        self.first_infeasible_position(ray)
            .unwrap_or_else(|| self.row_count() + 1)
    }

    fn sync_iteration_with_added(&mut self) {
        let iter = self
            .core
            .weakly_added_halfspaces
            .cardinality()
            .min(self.row_count());
        self.core.iter_state.iteration = iter;
    }

    #[cfg(debug_assertions)]
    fn debug_assert_iteration_prefix(&self) {
        let iter = self.core.iter_state.iteration.min(self.row_count());
        debug_assert_eq!(
            iter,
            self.core.weakly_added_halfspaces.cardinality(),
            "iteration {} diverges from weakly-added count {}",
            iter,
            self.core.weakly_added_halfspaces.cardinality()
        );
        for pos in 0..iter {
            let row = self.core.ctx.order_vector[pos];
            debug_assert!(
                self.core.weakly_added_halfspaces.contains(row),
                "row {} at position {} not marked weakly added",
                row,
                pos
            );
        }
        for row in 0..self.row_count() {
            if self.core.weakly_added_halfspaces.contains(row) {
                let in_prefix = self
                    .position_of_row(row)
                    .map(|pos| pos < iter)
                    .unwrap_or(false);
                debug_assert!(
                    in_prefix,
                    "weakly-added row {} not placed in iteration prefix",
                    row
                );
            }
        }
    }

    pub(crate) fn add_new_halfspace_dynamic(&mut self, row: Row) -> Result<(), Error>
    where
        ZR: crate::dd::mode::PreorderedBackend,
    {
        self.tableau_entering_for_row(row);
        let partition = self.evaluate_row_partition(row);
        let partition_view = Self::partition_view(&partition);

        if matches!(R::KIND, RepresentationKind::Inequality)
            && partition_view.negative.is_empty()
            && self.core.ctx.equality_kinds[row] != hullabaloo::types::InequalityKind::Equality
        {
            self.core.redundant_halfspaces.insert(row);
            self.core.weakly_added_halfspaces.insert(row);
            self.sync_iteration_with_added();
            self.sync_tableau_flag_for_row(row);
            if self.core.ray_graph.active_len() == self.core.ray_graph.weakly_feasible_len() {
                self.core.comp_status = ComputationStatus::AllFound;
            }
            self.recycle_partition(partition);
            return Ok(());
        }
        if partition_view.positive.is_empty() && partition_view.zero.is_empty() {
            self.core.ray_graph.deactivate_all();
            self.umpire.reset_infeasible_counts(self.row_count());
            self.clear_ray_indices();
            self.core.pending_new_rays.clear();
            self.core.comp_status = ComputationStatus::RegionEmpty;
            self.recycle_partition(partition);
            return Ok(());
        }

        if partition_view.negative.is_empty() {
            self.core.weakly_added_halfspaces.insert(row);
            self.sync_iteration_with_added();
            self.sync_tableau_flag_for_row(row);
            if self.core.ray_graph.active_len() == self.core.ray_graph.weakly_feasible_len() {
                self.core.comp_status = ComputationStatus::AllFound;
            }
            self.recycle_partition(partition);
            return Ok(());
        }

        let (new_id, inserted) = ZR::ensure_id_for_row(&mut self.core.ctx, row);
        if inserted {
            let use_incidence = !self.core.options.assumes_nondegeneracy()
                && ZR::USE_INCIDENCE_INDEX_FOR_CANDIDATE_TEST;
            let domain = ZR::domain_size(&self.core.ctx);
            for &ray_id in partition_view.zero {
                let Some(old_sig) = self.ray(ray_id).map(|ray_data| ray_data.zero_set_signature())
                else {
                    continue;
                };
                self.core.ray_index.unregister(ray_id, old_sig);
                let Some(new_sig) = self.ray_mut(ray_id).map(|ray_data| {
                    ray_data.zero_set_mut().insert(new_id);
                    ray_data.set_zero_set_count(ray_data.zero_set_count() + 1);
                    let sig = ray_data.zero_set().signature_u64();
                    ray_data.set_zero_set_signature(sig);
                    sig
                }) else {
                    continue;
                };
                self.core.ray_index.register(ray_id, new_sig);
                if use_incidence {
                    self.core
                        .ray_incidence
                        .add_zero_row_membership(ray_id, new_id, domain);
                }
            }
        }
        self.insert_added_halfspace_zero(row);

        self.with_active_ray_ids(|state, adjacency_candidates| {
            let use_bitset = ZR::USE_INCIDENCE_INDEX_FOR_CANDIDATE_TEST && {
                let ray_capacity = state.core.ray_incidence.ray_capacity();
                adjacency_candidates
                    .iter()
                    .all(|id| id.as_index() < ray_capacity)
            };

            if use_bitset {
                let ray_capacity = state.core.ray_incidence.ray_capacity();
                let mut candidate_set =
                    std::mem::replace(&mut state.core.candidate_ray_set, RowSet::new(0));
                if candidate_set.len() != ray_capacity {
                    candidate_set.resize(ray_capacity);
                }
                candidate_set.clear();
                for &rid in adjacency_candidates {
                    candidate_set.insert(rid.as_index());
                }
                for &pos in partition_view.positive {
                    for &neg in partition_view.negative {
                        let need = state.core.ray_incidence.ray_capacity();
                        if candidate_set.len() != need {
                            candidate_set.resize(need);
                        }
                        if state
                            .check_adjacency_common_with_candidates(
                                neg,
                                pos,
                                DegeneracyCandidates::BitSet(&candidate_set),
                            )
                            .is_none()
                        {
                            continue;
                        }
                        let _ = state.create_new_ray(neg, pos, row);
                    }
                }
                state.core.candidate_ray_set = candidate_set;
            } else {
                for &pos in partition_view.positive {
                    for &neg in partition_view.negative {
                        if state
                            .check_adjacency_common_with_candidates(
                                neg,
                                pos,
                                DegeneracyCandidates::Slice(adjacency_candidates),
                            )
                            .is_none()
                        {
                            continue;
                        }
                        let _ = state.create_new_ray(neg, pos, row);
                    }
                }
            }
        });
        self.delete_negative_rays(partition_view);
        self.core.added_halfspaces.insert(row);
        self.insert_added_halfspace_zero(row);
        self.prune_recent_rays();
        self.core.weakly_added_halfspaces.insert(row);
        self.sync_iteration_with_added();
        self.sync_tableau_flag_for_row(row);
        <UmpireHalfspaceMode<N, U, ZR> as HalfspaceMode>::on_halfspace_added(self, row);
        if self.core.ray_graph.active_len() == self.core.ray_graph.weakly_feasible_len() {
            self.core.comp_status = ComputationStatus::AllFound;
        }
        self.recycle_partition(partition);
        Ok(())
    }

    pub fn initialize_rays(&mut self)
    where
        ZR: crate::dd::mode::PreorderedBackend,
    {
        self.core.ray_graph.reset(self.row_count());
        self.umpire.reset_infeasible_counts(self.row_count());
        self.core.tableau.tableau_nonbasic = Self::init_nonbasic(self.col_count());
        self.core.tableau.tableau_basic_col_for_row = vec![-1; self.row_count()];
        self.add_artificial_ray();
        self.clear_ray_indices();

        self.core.initial_halfspaces.clear();
        for &row in &self.core.iter_state.initial_ray_index {
            if let Some(r) = row {
                self.core.initial_halfspaces.insert(r);
            }
        }

        self.core.added_halfspaces = self.core.initial_halfspaces.clone();
        self.core.weakly_added_halfspaces = self.core.initial_halfspaces.clone();
        for row_id in self.core.weakly_added_halfspaces.iter() {
            let row = row_id.as_index();
            let _ = ZR::ensure_id_for_row(&mut self.core.ctx, row);
        }
        self.rebuild_added_halfspaces_zero();
        self.sync_tableau_flags();
        let initial = self.core.initial_halfspaces.clone();
        self.umpire.bump_priority_rows(&mut self.core.ctx, &initial);
        self.core.edges_dirty = true;
        let last_row = self.core.ctx.order_vector.first().copied();
        let m = self.row_count();
        let wants_purify = self.umpire.wants_initial_purification();
        let mut expected_zero = RowSet::new(0);
        if wants_purify {
            expected_zero.resize(m);
        }
        for col in 0..self.col_count() {
            let mut vec = self
                .umpire
                .basis_column_vector(&self.core.tableau.basis, col);
            let pivot_row = self
                .core
                .iter_state
                .initial_ray_index
                .get(col)
                .copied()
                .flatten();
            if !self.umpire.normalize_vector(&mut vec) {
                continue;
            }
            if wants_purify {
                expected_zero.copy_from(&self.core.initial_halfspaces);
                expected_zero.union_inplace(&self.core.equality_set);
                if let Some(row) = pivot_row
                    && !self.core.equality_set.contains(row)
                {
                    expected_zero.remove(row);
                }

                if let Some(mut purified) =
                    self.umpire
                        .purify_vector_from_zero_set(&self.core.ctx, &expected_zero)
                    && self.umpire.normalize_vector(&mut purified)
                {
                    self.umpire.align_vector_in_place(&vec, &mut purified);
                    vec = purified;
                }
            }
            let has_initial = self
                .core
                .iter_state
                .initial_ray_index
                .get(col)
                .copied()
                .flatten()
                .is_some();
            let add_negative = !has_initial;
            let neg_vec = add_negative.then(|| {
                let mut neg = vec.clone();
                self.umpire.negate_vector_in_place(&mut neg);
                neg
            });
            // cddlib seeds use strict feasibility checks; do not apply relaxed handling here.
            let zero_set = self.take_zero_set();
            let mut ray_data = self
                .umpire
                .classify_vector(&self.core.ctx, vec, false, last_row, zero_set);
            let satisfies_equalities = self.core.equality_set.iter().all(|row_id| {
                let row = row_id.as_index();
                ZR::zero_set_contains_row(&self.core.ctx, ray_data.zero_set(), row)
            });
            if satisfies_equalities {
                self.umpire.on_ray_inserted(&self.core.ctx, &ray_data, false);
                let id = self.core.ray_graph.insert_active(ray_data);
                self.register_ray_id(id);
                if let Some(neg_vec) = neg_vec {
                    let zero_set = self.take_zero_set();
                    let neg_ray_data =
                        self.umpire
                            .classify_vector(&self.core.ctx, neg_vec, false, last_row, zero_set);
                    self.umpire
                        .on_ray_inserted(&self.core.ctx, &neg_ray_data, false);
                    let id = self.core.ray_graph.insert_active(neg_ray_data);
                    self.register_ray_id(id);
                }
            } else {
                self.recycle_zero_set_from_ray_data(&mut ray_data);
                self.umpire.recycle_ray_data(&mut ray_data);
            }
        }
        self.core.iter_state.iteration = self.col_count();
        self.with_active_ray_ids(|state, active_ids| {
            let mut floored = false;
            if state.core.iter_state.iteration < state.core.ctx.order_vector.len() {
                let floor_row = state.core.ctx.order_vector[state.core.iter_state.iteration];
                let floor_pos = state.core.iter_state.iteration;
                for &idx in active_ids {
                    if let Some(ray_pos) = state
                        .ray(idx)
                        .map(|ray| state.first_infeasible_position_or_m(ray))
                        && ray_pos < floor_pos
                        && let Some(ray_mut) = state.ray_mut(idx)
                    {
                        ray_mut.set_first_infeasible_row(Some(floor_row));
                        floored = true;
                    }
                }
            }
            if state.enforce_first_infeasible_floor(active_ids) {
                floored = true;
            }
            if floored {
                state.core.edges_dirty = true;
            }
        });
        <UmpireHalfspaceMode<N, U, ZR> as HalfspaceMode>::on_rays_initialized(self);
        self.core.iter_state.iteration = self.col_count().saturating_add(1);
        if self.core.iter_state.iteration > self.row_count() {
            self.core.iter_state.iteration = self.row_count();
        }
        if self.core.iter_state.iteration > self.row_count() {
            self.core.comp_status = ComputationStatus::AllFound;
        }

        // Start iteration counter aligned to rows that are actually marked added.
        // At initialization nothing has been added yet, so keep the iteration prefix empty.
        self.sync_iteration_with_added();
    }

    fn partition_view<'a>(partition: &'a RayPartitionOwned) -> RayPartition<'a> {
        RayPartition {
            negative: &partition.negative,
            positive: &partition.positive,
            zero: &partition.zero,
        }
    }

    pub(crate) fn recycle_partition(&mut self, mut partition: RayPartitionOwned) {
        partition.negative.clear();
        partition.positive.clear();
        partition.zero.clear();
        self.core.partitions = partition;
    }

    fn check_adjacency_common_with_candidates(
        &mut self,
        r1: RayId,
        r2: RayId,
        candidates: DegeneracyCandidates<'_>,
    ) -> Option<usize> {
        let required = self.adjacency_dimension().saturating_sub(2);
        let (ray1_zero, ray1_zero_count, ray2_zero, ray2_zero_count) = {
            let ray_graph = &self.core.ray_graph;
            let (Some(ray1), Some(ray2)) = (ray_graph.ray_data(r1), ray_graph.ray_data(r2))
            else {
                return None;
            };
            (
                ray1.zero_set(),
                ray1.zero_set_count(),
                ray2.zero_set(),
                ray2.zero_set_count(),
            )
        };

        let domain = ZR::domain_size(&self.core.ctx);
        let added_domain_full = self.core.added_halfspaces_zero_size == domain;
        let common = if added_domain_full {
            ray1_zero.count_intersection(ray2_zero)
        } else {
            self.core.adj_face.copy_from(ray1_zero);
            self.core
                .adj_face
                .intersection_two_inplace_and_count(ray2_zero, &self.core.added_halfspaces_zero)
        };

        if common < required {
            return None;
        }

        if self.core.options.assumes_nondegeneracy() {
            return Some(common);
        }

        let min_parent_saturators = if added_domain_full {
            ray1_zero_count.min(ray2_zero_count)
        } else {
            ray1_zero
                .count_intersection(&self.core.added_halfspaces_zero)
                .min(ray2_zero.count_intersection(&self.core.added_halfspaces_zero))
        };
        if min_parent_saturators == common.saturating_add(1) {
            return Some(common);
        }

        if added_domain_full {
            self.core.adj_face.copy_from(ray1_zero);
            self.core.adj_face.intersection_inplace(ray2_zero);
        }

        let contains = if ZR::USE_INCIDENCE_INDEX_FOR_CANDIDATE_TEST {
            match candidates {
                DegeneracyCandidates::Slice(indices) => self.core.ray_incidence.candidate_contains_face(
                    &self.core.adj_face,
                    domain,
                    indices,
                    r1,
                    r2,
                ),
                DegeneracyCandidates::BitSet(set) => self.core.ray_incidence.candidate_set_contains_face(
                    &self.core.adj_face,
                    domain,
                    set,
                    r1,
                    r2,
                ),
            }
        } else {
            self.candidate_contains_face_by_subset(&self.core.adj_face, common, candidates, r1, r2)
        };
        if contains {
            None
        } else {
            Some(common)
        }
    }

    #[inline]
    fn candidate_contains_face_by_subset(
        &self,
        face: &ZR::Set,
        face_cardinality: usize,
        candidates: DegeneracyCandidates<'_>,
        exclude_a: RayId,
        exclude_b: RayId,
    ) -> bool {
        match candidates {
            DegeneracyCandidates::Slice(indices) => {
                for &id in indices {
                    if id == exclude_a || id == exclude_b {
                        continue;
                    }
                    let Some(ray) = self.ray(id) else {
                        continue;
                    };
                    if ray.zero_set_count() < face_cardinality {
                        continue;
                    }
                    let candidate_zero = ray.zero_set();
                    if face.iter().all(|row_id| candidate_zero.contains(row_id)) {
                        return true;
                    }
                }
                false
            }
            DegeneracyCandidates::BitSet(set) => {
                for rid in set.iter() {
                    let id = RayId::from(rid.as_index());
                    if id == exclude_a || id == exclude_b {
                        continue;
                    }
                    let Some(ray) = self.ray(id) else {
                        continue;
                    };
                    if ray.zero_set_count() < face_cardinality {
                        continue;
                    }
                    let candidate_zero = ray.zero_set();
                    if face.iter().all(|row_id| candidate_zero.contains(row_id)) {
                        return true;
                    }
                }
                false
            }
        }
    }

    fn check_adjacency_common(
        &mut self,
        r1: RayId,
        r2: RayId,
        candidates: &[RayId],
    ) -> Option<usize> {
        self.check_adjacency_common_with_candidates(r1, r2, DegeneracyCandidates::Slice(candidates))
    }

    #[allow(dead_code)]
    pub(crate) fn check_adjacency(&mut self, r1: RayId, r2: RayId, candidates: &[RayId]) -> bool {
        self.check_adjacency_common(r1, r2, candidates).is_some()
    }

    #[allow(dead_code)]
    fn reclassify_active_rays(&mut self) {
        let mut ids = std::mem::take(&mut self.core.active_id_scratch);
        self.core.ray_graph.copy_active_ids(&mut ids);
        let relaxed = self.core.options.relaxed_enumeration();
        for id in ids.iter().copied() {
            let Some(ray_data) = self.core.ray_graph.ray_data_mut(id) else {
                continue;
            };
            self.umpire.reclassify_ray(&self.core.ctx, ray_data, relaxed);
        }
        self.core.active_id_scratch = ids;
        self.rebuild_ray_index();
        self.core.ray_graph.recompute_counts();
    }

    fn prune_recent_rays(&mut self) {
        self.core.pending_new_rays.clear();
    }

    pub(crate) fn delete_negative_rays(&mut self, partition: RayPartition<'_>) {
        let relaxed = self.core.options.relaxed_enumeration();
        for &idx in partition.negative.iter() {
            let ray_data = self.core.ray_graph.ray_data(idx);
            if let Some(ray_data) = ray_data {
                self.umpire
                    .on_ray_removed(&self.core.ctx, ray_data, relaxed);
            }
            self.unregister_ray_id(idx);
            if let Some(ray_data) = self.core.ray_graph.ray_data_mut(idx) {
                let mut set = ZR::empty_set(0);
                std::mem::swap(ray_data.zero_set_mut(), &mut set);
                self.core.zero_set_pool.push(set);
                self.umpire.recycle_ray_data(ray_data);
            }
        }
        self.core
            .ray_graph
            .remove_many_keep_order(partition.negative);
        self.discard_pending_rays(partition.negative);

        let mut zeros = std::mem::take(&mut self.core.active_id_scratch);
        zeros.clear();
        zeros.reserve(partition.zero.len() + self.core.pending_new_rays.len());
        zeros.extend_from_slice(partition.zero);

        let mut order = self.core.ray_graph.take_active_order();
        order.clear();
        order.reserve(partition.positive.len() + zeros.len() + self.core.pending_new_rays.len());
        order.extend_from_slice(partition.positive);

        let mut max_id = 0usize;
        for &rid in partition.positive.iter().chain(zeros.iter()) {
            max_id = max_id.max(rid.as_index());
        }
        for &rid in &self.core.pending_new_rays {
            max_id = max_id.max(rid.as_index());
        }
        if max_id >= self.core.removed_marks.len() {
            self.core.removed_marks.resize(max_id + 1, 0);
        }
        self.core.removed_epoch = self.core.removed_epoch.wrapping_add(1);
        if self.core.removed_epoch == 0 {
            self.core.removed_epoch = 1;
            self.core.removed_marks.fill(0);
        }
        let epoch = self.core.removed_epoch;
        for &rid in partition.positive.iter().chain(zeros.iter()) {
            self.core.removed_marks[rid.as_index()] = epoch;
        }

        for &id in &self.core.pending_new_rays {
            if self.core.removed_marks[id.as_index()] == epoch {
                continue;
            }
            let Some(ray) = self.ray(id) else { continue };
            if ray.last_sign() == Sign::Positive {
                order.push(id);
            } else {
                zeros.push(id);
            }
            self.core.removed_marks[id.as_index()] = epoch;
        }

        zeros.sort_unstable_by(|a, b| {
            let fa = self
                .ray(*a)
                .map(|r| self.first_infeasible_position_or_m(r))
                .unwrap_or(self.row_count());
            let fb = self
                .ray(*b)
                .map(|r| self.first_infeasible_position_or_m(r))
                .unwrap_or(self.row_count());
            fa.cmp(&fb).then_with(|| a.as_index().cmp(&b.as_index()))
        });

        let pos_head = order.first().copied();
        let pos_tail = order.last().copied();
        let zero_head = zeros.first().copied();
        let zero_tail = zeros.last().copied();
        let zero_count = zeros.len();
        order.extend_from_slice(&zeros);

        self.core
            .ray_graph
            .set_order_with_zero_count_unchecked(order, zero_count);

        self.core.lists = RayListHeads::default();
        self.core.lists.pos_head = pos_head;
        self.core.lists.pos_tail = pos_tail;
        self.core.lists.zero_head = zero_head;
        self.core.lists.zero_tail = zero_tail;
        self.core.active_id_scratch = zeros;
        if self.core.ray_graph.active_len() == 0 {
            self.core.comp_status = ComputationStatus::RegionEmpty;
            self.umpire.reset_infeasible_counts(self.row_count());
        }
    }

    pub(crate) fn refresh_edge_buckets(&mut self) {
        let row_cap = self.row_count();
        assert!(
            row_cap < usize::MAX,
            "row count overflow while refreshing edge buckets"
        );
        let target_len = self.core.ray_graph.edge_buckets_len().max(row_cap + 1);
        if target_len > 0 {
            self.core.ray_graph.ensure_edge_capacity(target_len - 1);
        }
        self.core.ray_graph.edge_count = 0;
        let buckets_to_process = std::mem::take(&mut self.core.ray_graph.non_empty_edge_buckets);
        for &bucket_idx in &buckets_to_process {
            debug_assert!(
                bucket_idx < self.core.ray_graph.edge_bucket_positions.len(),
                "edge bucket index out of range while clearing bucket positions"
            );
            self.core.ray_graph.edge_bucket_positions[bucket_idx] = None;
        }
        for bucket_idx in buckets_to_process {
            debug_assert!(
                bucket_idx < self.core.ray_graph.edges.len(),
                "edge bucket index out of range while draining buckets"
            );
            let mut bucket = Vec::new();
            std::mem::swap(&mut bucket, &mut self.core.ray_graph.edges[bucket_idx]);
            for edge in bucket.into_iter() {
                let target = self.edge_target_iteration(&edge);
                match target {
                    EdgeTarget::Scheduled(row) => self.core.ray_graph.schedule_edge(row, edge),
                    EdgeTarget::Stale(_) => {}
                    EdgeTarget::Discarded => {}
                }
            }
        }
    }

    pub(crate) fn edge_target_iteration(&self, edge: &AdjacencyEdge) -> EdgeTarget {
        let Some(ray1) = self.ray(edge.retained) else {
            return EdgeTarget::Discarded;
        };
        let Some(ray2) = self.ray(edge.removed) else {
            return EdgeTarget::Discarded;
        };
        let m = self.row_count();
        let f1 = self.first_infeasible_position_or_m(ray1);
        let f2 = self.first_infeasible_position_or_m(ray2);
        if f1 == f2 {
            return EdgeTarget::Discarded;
        }
        if f1 >= m && f2 >= m {
            return EdgeTarget::Discarded;
        }
        let fmin = f1.min(f2);
        if fmin < self.core.iter_state.iteration {
            return EdgeTarget::Stale(fmin);
        }
        EdgeTarget::Scheduled(fmin)
    }

    #[cfg(test)]
    pub(crate) fn add_ray(&mut self, vector: Vec<U::Scalar>) -> RayId {
        let relaxed = self.core.options.relaxed_enumeration();
        let last_row = self.core.ctx.order_vector.first().copied();
        let zero_set = self.take_zero_set();
        let ray_data = self
            .umpire
            .classify_vector(&self.core.ctx, vector, relaxed, last_row, zero_set);
        debug_assert!(!ray_data.is_feasible() || ray_data.is_weakly_feasible());
        self.umpire
            .on_ray_inserted(&self.core.ctx, &ray_data, relaxed);
        let id = self.core.ray_graph.insert_active(ray_data);
        self.register_ray_id(id);
        self.record_new_ray(id);
        id
    }

    pub(crate) fn add_artificial_ray(&mut self) -> RayId {
        if let Some(idx) = self.core.ray_graph.artificial_ray() {
            return idx;
        }
        let dimension = self.col_count().max(1);
        let vector = self.umpire.zero_vector(dimension);
        let last_row = self.core.ctx.order_vector.first().copied();
        let zero_set = self.take_zero_set();
        let ray_data = self
            .umpire
            .classify_vector(&self.core.ctx, vector, false, last_row, zero_set);
        let idx = self.core.ray_graph.insert_inactive(ray_data);
        self.core.ray_graph.set_artificial(idx);
        idx
    }

    pub(crate) fn evaluate_row_partition(&mut self, row: Row) -> RayPartitionOwned {
        self.core.lists = RayListHeads::default();
        let mut partition = std::mem::take(&mut self.core.partitions);
        partition.negative.clear();
        partition.positive.clear();
        partition.zero.clear();
        let mut active_order = self.core.ray_graph.take_active_order();
        {
            let ctx = &self.core.ctx;
            let (umpire, ray_graph) = (&mut self.umpire, &mut self.core.ray_graph);
            for &idx in active_order.iter() {
                let Some(ray_data) = ray_graph.ray_data_mut(idx) else {
                    continue;
                };
                let sign = umpire.classify_ray(ctx, ray_data, row);
                match sign {
                    Sign::Negative => partition.negative.push(idx),
                    Sign::Positive => partition.positive.push(idx),
                    Sign::Zero => partition.zero.push(idx),
                }
            }
        }

        self.core.lists.neg_head = partition.negative.first().copied();
        self.core.lists.neg_tail = partition.negative.last().copied();
        self.core.lists.pos_head = partition.positive.first().copied();
        self.core.lists.pos_tail = partition.positive.last().copied();
        self.core.lists.zero_head = partition.zero.first().copied();
        self.core.lists.zero_tail = partition.zero.last().copied();

        let floored_pos = self.enforce_first_infeasible_floor(&partition.positive);
        let floored_zero = self.enforce_first_infeasible_floor(&partition.zero);
        if floored_pos || floored_zero {
            self.core.edges_dirty = true;
        }

        let mut cursor = 0usize;
        for &id in partition
            .negative
            .iter()
            .chain(partition.positive.iter())
            .chain(partition.zero.iter())
        {
            if cursor < active_order.len() {
                active_order[cursor] = id;
            } else {
                active_order.push(id);
            }
            cursor += 1;
        }
        active_order.truncate(cursor);

        self.core
            .ray_graph
            .set_order_with_zero_count_unchecked(active_order, partition.zero.len());
        partition
    }

    pub(crate) fn create_new_ray(
        &mut self,
        r1: RayId,
        r2: RayId,
        row: Row,
    ) -> Option<RayId> {
        if r1 == r2 {
            return None;
        }

        let relaxed = self.core.options.relaxed_enumeration();
        let parent_a = self.core.ray_graph.ray_key(r1);
        let parent_b = self.core.ray_graph.ray_key(r2);
        let ray_data_res = {
            let zero_set = self.take_zero_set();
            if self.core.ray_graph.ray_data(r1).is_none() || self.core.ray_graph.ray_data(r2).is_none() {
                self.recycle_zero_set(zero_set);
                return None;
            }
            let ray1 = self.core.ray_graph.ray_data(r1).expect("ray id checked");
            let ray2 = self.core.ray_graph.ray_data(r2).expect("ray id checked");
            self.umpire.generate_new_ray(
                &self.core.ctx,
                (r1, ray1, r2, ray2),
                row,
                relaxed,
                zero_set,
            )
        };
        let mut ray_data = match ray_data_res {
            Ok(ray_data) => ray_data,
            Err(zero_set) => {
                self.recycle_zero_set(zero_set);
                return None;
            }
        };
        if self.ray_exists(&ray_data) {
            self.core.dedup_drops = self.core.dedup_drops.wrapping_add(1);
            self.recycle_zero_set_from_ray_data(&mut ray_data);
            self.umpire.recycle_ray_data(&mut ray_data);
            return None;
        }

        self.umpire.on_ray_inserted(&self.core.ctx, &ray_data, relaxed);
        let id = self.core.ray_graph.insert_active(ray_data);
        self.core
            .ray_graph
            .set_ray_origin(id, Some(parent_a), Some(parent_b), Some(row));
        self.register_ray_id(id);
        self.record_new_ray(id);
        Some(id)
    }

    pub fn expand_ray_vector(&self, compact: &[N]) -> Vec<N> {
        assert_eq!(
            compact.len(),
            self.col_count(),
            "expand_ray_vector expects a vector of length col_count"
        );
        assert_eq!(
            self.core.iter_state.newcol.len(),
            self.core.iter_state.d_orig,
            "column remapping must match original dimension"
        );
        if !self.core.iter_state.col_reduced || self.core.iter_state.d_orig == self.col_count() {
            return compact.to_vec();
        }
        let mut full = vec![N::zero(); self.core.iter_state.d_orig];
        for (orig_idx, &map) in self.core.iter_state.newcol.iter().enumerate() {
            if let Some(m) = map {
                assert!(m < compact.len(), "column remapping out of range");
                // SAFETY: `orig_idx < d_orig == full.len()` by construction, and `m < compact.len()` asserted above.
                unsafe {
                    *full.get_unchecked_mut(orig_idx) = compact.get_unchecked(m).clone();
                }
            }
        }
        full
    }

    pub(crate) fn record_new_ray(&mut self, id: RayId) {
        self.core.pending_new_rays.push(id);
    }

    fn discard_pending_rays(&mut self, removed: &[RayId]) {
        if self.core.pending_new_rays.is_empty() || removed.is_empty() {
            return;
        }
        self.core.removed_epoch = self.core.removed_epoch.wrapping_add(1);
        if self.core.removed_epoch == 0 {
            self.core.removed_epoch = 1;
            self.core.removed_marks.fill(0);
        }

        let mut max_id = 0usize;
        for &rid in removed {
            max_id = max_id.max(rid.as_index());
        }
        for &rid in &self.core.pending_new_rays {
            max_id = max_id.max(rid.as_index());
        }
        if max_id >= self.core.removed_marks.len() {
            self.core.removed_marks.resize(max_id + 1, 0);
        }

        let epoch = self.core.removed_epoch;
        for &rid in removed {
            self.core.removed_marks[rid.as_index()] = epoch;
        }

        self.core
            .pending_new_rays
            .retain(|rid| self.core.removed_marks[rid.as_index()] != epoch);
    }

    pub(crate) fn ray(&self, index: RayId) -> Option<&U::RayData> {
        self.core.ray_graph.ray_data(index)
    }

    pub(crate) fn ray_mut(&mut self, index: RayId) -> Option<&mut U::RayData> {
        self.core.ray_graph.ray_data_mut(index)
    }

    pub(crate) fn with_active_ray_ids<T>(&mut self, f: impl FnOnce(&mut Self, &[RayId]) -> T) -> T {
        let mut ids = std::mem::take(&mut self.core.active_id_scratch);
        self.core.ray_graph.copy_active_ids(&mut ids);
        let result = f(self, &ids);
        self.core.active_id_scratch = ids;
        result
    }

    pub(crate) fn sync_tableau_flags(&mut self) {
        let m = self.row_count();
        if self.core.tableau.row_status.len() != m {
            self.core.tableau.row_status = vec![-1; m];
        } else {
            self.core.tableau.row_status.fill(-1);
        }
        for row in self.core.weakly_added_halfspaces.iter() {
            self.core.tableau.row_status[row.as_index()] = 0;
        }
        for row in self.core.added_halfspaces.iter() {
            self.core.tableau.row_status[row.as_index()] = 1;
        }
    }

    #[inline]
    pub(crate) fn sync_tableau_flag_for_row(&mut self, row: Row) {
        if row >= self.row_count() {
            return;
        }
        if self.core.tableau.row_status.len() != self.row_count() {
            self.sync_tableau_flags();
            return;
        }
        let flag = if self.core.added_halfspaces.contains(row) {
            1
        } else if self.core.weakly_added_halfspaces.contains(row) {
            0
        } else {
            -1
        };
        self.core.tableau.row_status[row] = flag;
    }

    pub(crate) fn enforce_first_infeasible_floor(&mut self, indices: &[RayId]) -> bool {
        let min_order = self.core.iter_state.iteration;
        if min_order >= self.core.ctx.order_vector.len() {
            return false;
        }
        let floor_row = self.core.ctx.order_vector[min_order];
        let mut changed = false;
        for &idx in indices {
            let should_floor = match self.ray(idx) {
                Some(ray) => {
                    if ray.last_sign() == Sign::Negative {
                        false
                    } else {
                        self.first_infeasible_position(ray)
                            .map(|pos| pos < min_order)
                            .unwrap_or(false)
                    }
                }
                None => false,
            };
            if should_floor && let Some(ray) = self.ray_mut(idx) {
                ray.set_first_infeasible_row(Some(floor_row));
                changed = true;
            }
        }
        changed
    }

    fn adjacency_dimension(&self) -> usize {
        self.col_count()
    }
}

impl<N: Num, R: Representation, U> ConeEngine<N, R, U, crate::dd::RowRepr>
where
    U: Umpire<N, crate::dd::RowRepr>,
{
    pub(crate) fn add_new_halfspace_preordered(&mut self, row: Row) -> Result<(), Error> {
        let iter_pos = self
            .order_position(row)
            .unwrap_or_else(|| panic!("row {row} not present in order vector"));
        self.core.iter_state.iteration = iter_pos;
        self.tableau_entering_for_row(row);

        let partition = self.evaluate_row_partition(row);
        let partition_view = Self::partition_view(&partition);

        if matches!(R::KIND, RepresentationKind::Inequality)
            && partition_view.negative.is_empty()
            && self.core.ctx.equality_kinds[row] != hullabaloo::types::InequalityKind::Equality
        {
            self.core.redundant_halfspaces.insert(row);
            self.core.weakly_added_halfspaces.insert(row);
            self.sync_iteration_with_added();
            self.sync_tableau_flag_for_row(row);
            if self.core.ray_graph.active_len() == self.core.ray_graph.weakly_feasible_len() {
                self.core.comp_status = ComputationStatus::AllFound;
            }
            self.recycle_partition(partition);
            return Ok(());
        }
        if partition_view.positive.is_empty() && partition_view.zero.is_empty() {
            self.core.ray_graph.deactivate_all();
            self.umpire.reset_infeasible_counts(self.row_count());
            self.clear_ray_indices();
            self.core.pending_new_rays.clear();
            self.core.comp_status = ComputationStatus::RegionEmpty;
            self.recycle_partition(partition);
            return Ok(());
        }

        let valid_first = Some(partition_view.positive);
        self.process_iteration_edges(iter_pos, self.core.ctx.order_vector[iter_pos], valid_first);
        self.delete_negative_rays(partition_view);
        self.prune_recent_rays();
        self.core.added_halfspaces.insert(row);
        self.insert_added_halfspace_zero(row);
        self.core.weakly_added_halfspaces.insert(row);
        self.sync_iteration_with_added();
        self.sync_tableau_flag_for_row(row);

        if self.core.iter_state.iteration < self.row_count() && self.core.ray_graph.zero_len() > 1 {
            let zero_rays = std::mem::take(&mut self.core.active_id_scratch);
            if zero_rays.len() > 1 {
                self.update_edges(&zero_rays);
            }
            self.core.active_id_scratch = zero_rays;
        }

        <UmpireHalfspaceMode<N, U, crate::dd::RowRepr> as HalfspaceMode>::on_halfspace_added(
            self, row,
        );
        if self.core.ray_graph.active_len() == self.core.ray_graph.weakly_feasible_len() {
            self.core.comp_status = ComputationStatus::AllFound;
        }
        self.recycle_partition(partition);
        Ok(())
    }

    fn conditional_add_edge(&mut self, r1: RayId, r2: RayId, candidates: DegeneracyCandidates<'_>) {
        let f1 = match self.ray(r1) {
            Some(ray) => self.first_infeasible_position_or_m(ray),
            None => return,
        };
        let f2 = match self.ray(r2) {
            Some(ray) => self.first_infeasible_position_or_m(ray),
            None => return,
        };
        let (fmin, rmin, rmax) = if f1 <= f2 { (f1, r1, r2) } else { (f2, r2, r1) };
        if fmin >= self.core.ctx.order_vector.len() {
            return;
        }
        if self.ray(rmax).is_none() {
            return;
        }

        let max_zero = {
            let ray_graph = &self.core.ray_graph;
            let Some(ray) = ray_graph.ray_data(rmax) else {
                return;
            };
            ray.zero_set()
        };
        self.core.adj_face.copy_from(max_zero);

        let fmin_row = self
            .core
            .ctx
            .order_vector
            .get(fmin)
            .copied()
            .unwrap_or(self.row_count());
        let strong_sep = !self.core.adj_face.contains(fmin_row);
        if !strong_sep {
            return;
        }

        let min_zero = {
            let ray_graph = &self.core.ray_graph;
            let Some(ray) = ray_graph.ray_data(rmin) else {
                return;
            };
            ray.zero_set()
        };
        self.core.adj_face.intersection_inplace(min_zero);

        self.core.count_intersections += 1;
        let mut last_chance = true;
        for iteration in (self.core.iter_state.iteration + 1)..fmin {
            let row_idx = self.core.ctx.order_vector[iteration];
            let contains_row = self.core.adj_face.contains(row_idx);
            if contains_row && !self.core._strict_inequality_set.contains(row_idx) {
                last_chance = false;
                self.core.count_intersections_bad += 1;
                break;
            }
        }
        if last_chance {
            self.core.count_intersections_good += 1;
            let use_added = !self.core.added_halfspaces.is_empty();
            if use_added {
                self.core
                    .adj_face
                    .intersection_inplace(&self.core.added_halfspaces);
            } else {
                self.core.adj_face.intersection_inplace(&self.core.ground_set);
            }
            let required = self.adjacency_dimension().saturating_sub(2);
            let mut adjacent = self.core.adj_face.cardinality_at_least(required);
            if adjacent && !self.core.options.assumes_nondegeneracy() {
                let domain = crate::dd::RowRepr::domain_size(&self.core.ctx);
                let contains = match candidates {
                    DegeneracyCandidates::Slice(indices) => self.core.ray_incidence.candidate_contains_face(
                        &self.core.adj_face,
                        domain,
                        indices,
                        r1,
                        r2,
                    ),
                    DegeneracyCandidates::BitSet(set) => self.core.ray_incidence.candidate_set_contains_face(
                        &self.core.adj_face,
                        domain,
                        set,
                        r1,
                        r2,
                    ),
                };
                if contains {
                    adjacent = false;
                }
            }
            if adjacent {
                let edge = AdjacencyEdge {
                    retained: rmax,
                    removed: rmin,
                };
                self.core.ray_graph.queue_edge(fmin, edge);
            }
        }
    }

    pub(crate) fn create_initial_edges(&mut self) {
        self.core.iter_state.iteration = self.col_count();
        self.with_active_ray_ids(|state, rays| {
            if rays.len() < 2 {
                return;
            }
            for i in 0..rays.len() - 1 {
                let r1 = rays[i];
                let f1 = state
                    .ray(r1)
                    .map(|r| state.first_infeasible_position_or_m(r))
                    .unwrap_or(state.row_count());
                for &r2 in rays.iter().skip(i + 1) {
                    let f2 = state
                        .ray(r2)
                        .map(|r| state.first_infeasible_position_or_m(r))
                        .unwrap_or(state.row_count());
                    if f1 == f2 {
                        continue;
                    }
                    state.conditional_add_edge(r1, r2, DegeneracyCandidates::Slice(rays));
                }
            }
        });
    }

    pub(crate) fn update_edges(&mut self, zero_rays: &[RayId]) {
        if zero_rays.len() < 2 {
            return;
        }
        for i in 0..zero_rays.len() - 1 {
            let r1 = zero_rays[i];
            let fi = self
                .ray(r1)
                .map(|r| self.first_infeasible_position_or_m(r))
                .unwrap_or(self.row_count());
            for &r2 in zero_rays.iter().skip(i + 1) {
                let f2 = self
                    .ray(r2)
                    .map(|r| self.first_infeasible_position_or_m(r))
                    .unwrap_or(self.row_count());
                if f2 <= fi {
                    continue;
                }
                self.conditional_add_edge(r1, r2, DegeneracyCandidates::Slice(zero_rays));
            }
        }
    }

    fn process_iteration_edges(&mut self, iteration: Row, row: Row, valid_first: Option<&[RayId]>) {
        if iteration >= self.core.ray_graph.edge_buckets_len() {
            return;
        }
        if self.core.edges_dirty {
            self.refresh_edge_buckets();
            self.core.edges_dirty = false;
        }

        let valid_first = valid_first.unwrap_or(&[]);
        let edges = self.core.ray_graph.take_edges(iteration);

        if self.core.options.assumes_nondegeneracy() {
            for edge in edges {
                let f1 = match self.ray(edge.retained) {
                    Some(ray) => self.first_infeasible_position_or_m(ray),
                    None => continue,
                };

                let Some(new_idx) = self.create_new_ray(edge.removed, edge.retained, row) else {
                    continue;
                };

                let f2 = match self.ray(new_idx) {
                    Some(ray) => self.first_infeasible_position_or_m(ray),
                    None => continue,
                };

                if f1 != f2 {
                    self.conditional_add_edge(
                        edge.retained,
                        new_idx,
                        DegeneracyCandidates::Slice(valid_first),
                    );
                }
            }
            return;
        }

        let mut candidate_set = std::mem::replace(&mut self.core.candidate_ray_set, RowSet::new(0));
        candidate_set.resize(self.core.ray_incidence.ray_capacity());
        candidate_set.clear();
        for &rid in valid_first {
            candidate_set.insert(rid.as_index());
        }

        for edge in edges {
            let f1 = match self.ray(edge.retained) {
                Some(ray) => self.first_infeasible_position_or_m(ray),
                None => continue,
            };

            let Some(new_idx) = self.create_new_ray(edge.removed, edge.retained, row) else {
                continue;
            };

            let f2 = match self.ray(new_idx) {
                Some(ray) => self.first_infeasible_position_or_m(ray),
                None => continue,
            };

            if f1 != f2 {
                let need = self.core.ray_incidence.ray_capacity();
                if candidate_set.len() != need {
                    candidate_set.resize(need);
                }
                self.conditional_add_edge(
                    edge.retained,
                    new_idx,
                    DegeneracyCandidates::BitSet(&candidate_set),
                );
            }
        }
        self.core.candidate_ray_set = candidate_set;
    }
}

impl<N: Num, R: Representation, U, ZR> ConeEngine<N, R, U, ZR>
where
    U: Umpire<N, ZR>,
    ZR: ZeroRepr,
{
    pub fn feasibility_indices(&mut self, row: Row) -> (usize, usize) {
        let active = self.core.ray_graph.active_len();
        if row >= self.row_count() {
            return (active, 0);
        }

        let mut infeasible = 0usize;
        let mut ids = std::mem::take(&mut self.core.active_id_scratch);
        self.core.ray_graph.copy_active_ids(&mut ids);
        for id in ids.iter().copied() {
            let Some(ray_data) = self.core.ray_graph.ray_data(id) else {
                continue;
            };
            if self
                .umpire
                .sign_for_row_on_ray(&self.core.ctx, ray_data, row)
                == Sign::Negative
            {
                infeasible += 1;
            }
        }
        self.core.active_id_scratch = ids;

        assert!(active >= infeasible, "infeasible count exceeds active rays");
        (active - infeasible, infeasible)
    }
}
