use bone_solver::{
    ConstraintSystem, LineHandle, NewtonConfig, PointHandle, Residual, SolverError,
    assemble_jacobian, decompose, evaluate_residuals, jacobian_triplets, residual_norm,
    solve_newton, solve_newton_decomposed, sparsity_pattern,
};
use bone_types::{NewtonDamping, Parameter, ParameterIndex, ResidualIndex, SolverResidual};

fn p(idx: u32) -> ParameterIndex {
    ParameterIndex::new(idx)
}

fn point(x: u32, y: u32) -> PointHandle {
    PointHandle { x: p(x), y: p(y) }
}

fn parameters(values: &[f64]) -> Vec<Parameter> {
    values.iter().copied().map(Parameter::new).collect()
}

#[test]
fn horizontal_residual_row_and_jacobian_match_hand_derivation() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 1.5, 3.0, 4.0]),
        vec![Residual::Horizontal(LineHandle {
            a: point(0, 1),
            b: point(2, 3),
        })],
    );
    assert_eq!(system.row_count(), 1);
    let values = evaluate_residuals(&system, &[0.0, 1.5, 3.0, 4.0]);
    assert!((values[0] - (1.5 - 4.0)).abs() < 1e-15);
    let triplets = jacobian_triplets(&system, &[0.0, 1.5, 3.0, 4.0]);
    assert_eq!(triplets.len(), 2);
    let zero = ResidualIndex::new(0);
    assert!(triplets.contains(&(zero, p(1), 1.0)));
    assert!(triplets.contains(&(zero, p(3), -1.0)));
}

#[test]
fn pin_solves_in_one_newton_step() {
    let system = ConstraintSystem::new(
        parameters(&[3.0]),
        vec![Residual::Pin {
            param: p(0),
            target: 7.0,
        }],
    );
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("pin is linear, must converge");
    };
    assert!((out[0].value() - 7.0).abs() < 1e-9);
}

#[test]
fn coincident_point_point_is_two_row_residual() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0, 1.0, 2.0]),
        vec![Residual::CoincidentPointPoint(point(0, 1), point(2, 3))],
    );
    assert_eq!(system.row_count(), 2);
    let r = evaluate_residuals(&system, &[0.0, 0.0, 1.0, 2.0]);
    assert_eq!(r.len(), 2);
    assert!((r[0] - (0.0 - 1.0)).abs() < 1e-15);
    assert!((r[1] - (0.0 - 2.0)).abs() < 1e-15);
}

#[test]
fn sparsity_pattern_drops_duplicate_positions() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0]),
        vec![Residual::Pin {
            param: p(0),
            target: 1.0,
        }],
    );
    let pattern = sparsity_pattern(&system);
    assert_eq!(pattern.rows(), 1);
    assert_eq!(pattern.cols(), 2);
    assert_eq!(pattern.entries().len(), 1);
}

#[test]
fn default_damping_regularises_rank_deficient_jacobian() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 3.0, 5.0, 7.0]),
        vec![Residual::Horizontal(LineHandle {
            a: point(0, 1),
            b: point(2, 3),
        })],
    );
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("default damping must regularise rank-deficient J^T J");
    };
    assert!((out[1].value() - out[3].value()).abs() < 1e-9);
}

#[test]
fn zero_damping_on_rank_deficient_jacobian_does_not_silently_succeed() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 3.0, 5.0, 7.0]),
        vec![Residual::Horizontal(LineHandle {
            a: point(0, 1),
            b: point(2, 3),
        })],
    );
    let cfg = NewtonConfig {
        damping: NewtonDamping::new(0.0),
        ..NewtonConfig::DEFAULT
    };
    let Err(err) = solve_newton(&system, cfg) else {
        panic!("damping=0 on rank-deficient system must surface an error");
    };
    assert!(
        matches!(
            err,
            SolverError::InvalidSolutionFound { .. } | SolverError::NoSolutionFound { .. }
        ),
        "expected InvalidSolutionFound or NoSolutionFound, got {err:?}"
    );
}

#[test]
fn residual_norm_matches_euclidean_length() {
    let out = evaluate_residuals(
        &ConstraintSystem::new(
            parameters(&[3.0]),
            vec![Residual::Pin {
                param: p(0),
                target: 0.0,
            }],
        ),
        &[3.0],
    );
    let norm = residual_norm(&out);
    assert!((norm.value() - 3.0).abs() < 1e-15);
    let _ = SolverResidual::new(norm.value());
}

#[test]
fn decomposed_singular_index_is_translated_to_global_parameter() {
    let pins: Vec<Residual> = (0..10)
        .map(|i| Residual::Pin {
            param: p(i),
            target: 0.0,
        })
        .collect();
    let residuals: Vec<Residual> = pins
        .into_iter()
        .chain(std::iter::once(Residual::Horizontal(LineHandle {
            a: point(10, 11),
            b: point(12, 13),
        })))
        .collect();
    let seed = [0.0; 10]
        .iter()
        .copied()
        .chain([0.0, 5.0, 1.0, 7.0])
        .collect::<Vec<_>>();
    let system = ConstraintSystem::new(parameters(&seed), residuals);
    let decomp = decompose(&system);
    let cfg = NewtonConfig {
        damping: NewtonDamping::new(0.0),
        ..NewtonConfig::DEFAULT
    };
    match solve_newton_decomposed(&system, &decomp, cfg) {
        Err(SolverError::InvalidSolutionFound { at }) => {
            assert!(
                at.as_usize() >= 10,
                "expected global parameter index (>=10), got {at:?}"
            );
        }
        Err(SolverError::NoSolutionFound { .. }) => {}
        other => panic!("expected InvalidSolutionFound or NoSolutionFound, got {other:?}"),
    }
}

#[test]
fn symmetric_residual_is_zero_when_b_is_reflection_of_a() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 1.0, 0.0, -1.0, -2.0, 0.0, 2.0, 0.0]),
        vec![Residual::Symmetric {
            a: point(0, 1),
            b: point(2, 3),
            axis: LineHandle {
                a: point(4, 5),
                b: point(6, 7),
            },
        }],
    );
    assert_eq!(system.row_count(), 2);
    let values = evaluate_residuals(&system, &[0.0, 1.0, 0.0, -1.0, -2.0, 0.0, 2.0, 0.0]);
    assert!(values[0].abs() < 1e-15, "perpendicularity: {}", values[0]);
    assert!(values[1].abs() < 1e-15, "midpoint-on-axis: {}", values[1]);
}

#[test]
fn symmetric_residual_is_nonzero_when_b_drifts_off_reflection() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 1.0, 0.5, -1.0, -2.0, 0.0, 2.0, 0.0]),
        vec![Residual::Symmetric {
            a: point(0, 1),
            b: point(2, 3),
            axis: LineHandle {
                a: point(4, 5),
                b: point(6, 7),
            },
        }],
    );
    let values = evaluate_residuals(&system, &[0.0, 1.0, 0.5, -1.0, -2.0, 0.0, 2.0, 0.0]);
    assert!(
        values[0].abs() > 0.0 || values[1].abs() > 0.0,
        "drift must surface in at least one residual",
    );
}

#[test]
fn symmetric_solves_to_b_as_reflection_of_a_about_x_axis() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 1.0, 0.5, 1.5, -2.0, 0.0, 2.0, 0.0]),
        vec![
            Residual::Pin {
                param: p(0),
                target: 0.0,
            },
            Residual::Pin {
                param: p(1),
                target: 1.0,
            },
            Residual::Pin {
                param: p(4),
                target: -2.0,
            },
            Residual::Pin {
                param: p(5),
                target: 0.0,
            },
            Residual::Pin {
                param: p(6),
                target: 2.0,
            },
            Residual::Pin {
                param: p(7),
                target: 0.0,
            },
            Residual::Symmetric {
                a: point(0, 1),
                b: point(2, 3),
                axis: LineHandle {
                    a: point(4, 5),
                    b: point(6, 7),
                },
            },
        ],
    );
    let Ok(out) = solve_newton(&system, NewtonConfig::DEFAULT) else {
        panic!("symmetric over fixed axis must converge");
    };
    assert!(
        (out[2].value() - 0.0).abs() < 1e-6,
        "b.x = a.x for reflection over x-axis: got {}",
        out[2].value()
    );
    assert!(
        (out[3].value() - -1.0).abs() < 1e-6,
        "b.y = -a.y for reflection over x-axis: got {}",
        out[3].value()
    );
}

#[test]
fn assemble_jacobian_and_triplets_agree_on_dense_shape() {
    let system = ConstraintSystem::new(
        parameters(&[0.0, 0.0, 2.0, 0.0]),
        vec![Residual::Horizontal(LineHandle {
            a: point(0, 1),
            b: point(2, 3),
        })],
    );
    let params = vec![0.0, 0.0, 2.0, 0.0];
    let dense = assemble_jacobian(&system, &params);
    assert_eq!(dense.nrows(), 1);
    assert_eq!(dense.ncols(), 4);
    let triplets = jacobian_triplets(&system, &params);
    let sum: f64 = triplets
        .iter()
        .map(|(r, c, v)| dense[(r.as_usize(), c.as_usize())] - v)
        .sum();
    assert!(sum.abs() < 1e-15);
}
