#![cfg(feature = "rug")]

use anyhow::Result;
use calculo::num::Num;
use cddlib_rs::{Matrix as CddMatrix, NumberType, Polyhedron as CddPolyhedron};
use howzat::lrs::Options;
use howzat::matrix::LpMatrixBuilder;
use howzat::polyhedron::{LrsConfig, PolyhedronOptions, PolyhedronOutput};
use hullabaloo::types::{Generator, IncidenceOutput};

#[derive(Clone, Debug)]
enum GeneratorRow {
    Vertex(Vec<f64>),
    Ray(Vec<f64>),
}

impl GeneratorRow {
    fn is_vertex(&self) -> bool {
        matches!(self, GeneratorRow::Vertex(_))
    }

    fn coords(&self) -> &[f64] {
        match self {
            GeneratorRow::Vertex(coords) | GeneratorRow::Ray(coords) => coords,
        }
    }
}

fn build_generator_matrix(
    rows: &[GeneratorRow],
) -> Result<howzat::matrix::LpMatrix<f64, Generator>> {
    let rows: Vec<Vec<f64>> = rows
        .iter()
        .map(|g| {
            let coords = g.coords();
            let mut row = Vec::with_capacity(coords.len() + 1);
            row.push(if g.is_vertex() { 1.0 } else { 0.0 });
            row.extend_from_slice(coords);
            row
        })
        .collect();
    Ok(LpMatrixBuilder::from_rows(rows).build())
}

fn cdd_incidence(rows: &[GeneratorRow]) -> Result<Vec<Vec<usize>>> {
    let mut cdd: CddMatrix<f64, Generator> = CddMatrix::new(
        rows.len(),
        rows.first()
            .map(|r| r.coords().len() + 1)
            .unwrap_or_default(),
        NumberType::Real,
    )?;
    for (r, g) in rows.iter().enumerate() {
        let coords = g.coords();
        cdd.set_real(r, 0, if g.is_vertex() { 1.0 } else { 0.0 });
        for (c, &val) in coords.iter().enumerate() {
            cdd.set_real(r, c + 1, val);
        }
        cdd.set_generator_type(r, g.is_vertex());
    }
    let poly = CddPolyhedron::from_matrix(&cdd)?;
    let mut facets = poly.incidence()?.to_adjacency_lists();
    for face in facets.iter_mut() {
        face.sort_unstable();
    }
    facets.sort();
    Ok(facets)
}

fn howzat_incidence_lrs(rows: &[GeneratorRow]) -> Result<Vec<Vec<usize>>> {
    let matrix = build_generator_matrix(rows)?;
    let input_row_count = matrix.row_count();
    let eps = f64::default_eps();
    let poly =
        PolyhedronOutput::<f64, Generator>::from_matrix_lrs_as_exact::<calculo::num::RugRat, f64>(
            matrix,
            LrsConfig {
                poly: PolyhedronOptions {
                    output_incidence: IncidenceOutput::Set,
                    ..PolyhedronOptions::default()
                },
                lrs: Options::default(),
            },
            &eps,
        )?;
    let incidence = poly.incidence_required();
    let linearity = poly.output().linearity().clone();
    let mut faces = Vec::new();
    for idx in 0..incidence.family_size() {
        if linearity.contains(idx) {
            continue;
        }
        let set = incidence
            .set(idx)
            .cloned()
            .unwrap_or_else(|| hullabaloo::types::RowSet::new(input_row_count));
        let mut verts: Vec<usize> = set.iter().map(|v| v.as_index()).collect();
        verts.sort_unstable();
        faces.push(verts);
    }
    faces.sort();
    Ok(faces)
}

fn verify_against_cdd(rows: &[GeneratorRow]) -> Result<()> {
    let expected = cdd_incidence(rows)?;
    let actual = howzat_incidence_lrs(rows)?;
    assert_eq!(
        actual, expected,
        "howzat(LRS) incidence differed from cddlib for generators {:?}",
        rows
    );
    Ok(())
}

#[test]
fn triangle_matches_cdd_with_lrs() -> Result<()> {
    let rows = vec![
        GeneratorRow::Vertex(vec![0.0, 0.0]),
        GeneratorRow::Vertex(vec![1.0, 0.0]),
        GeneratorRow::Vertex(vec![0.0, 1.0]),
    ];
    verify_against_cdd(&rows)
}

#[test]
fn triangular_prism_matches_cdd_with_lrs() -> Result<()> {
    let rows = vec![
        GeneratorRow::Vertex(vec![0.0, 0.0, 0.0]),
        GeneratorRow::Vertex(vec![1.0, 0.0, 0.0]),
        GeneratorRow::Vertex(vec![0.0, 1.0, 0.0]),
        GeneratorRow::Vertex(vec![0.0, 0.0, 1.0]),
        GeneratorRow::Vertex(vec![1.0, 0.0, 1.0]),
        GeneratorRow::Vertex(vec![0.0, 1.0, 1.0]),
    ];
    verify_against_cdd(&rows)
}

#[test]
fn cube_matches_cdd_with_lrs() -> Result<()> {
    let mut rows = Vec::new();
    for &x in &[0.0, 1.0] {
        for &y in &[0.0, 1.0] {
            for &z in &[0.0, 1.0] {
                rows.push(GeneratorRow::Vertex(vec![x, y, z]));
            }
        }
    }
    verify_against_cdd(&rows)
}

#[test]
fn positive_orthant_cone_matches_cdd_with_lrs() -> Result<()> {
    // Cone generated by the standard basis in R^3, with an explicit apex vertex at the origin.
    // Facets are x>=0, y>=0, z>=0; no lineality.
    let rows = vec![
        GeneratorRow::Vertex(vec![0.0, 0.0, 0.0]),
        GeneratorRow::Ray(vec![1.0, 0.0, 0.0]),
        GeneratorRow::Ray(vec![0.0, 1.0, 0.0]),
        GeneratorRow::Ray(vec![0.0, 0.0, 1.0]),
    ];
    verify_against_cdd(&rows)
}
