// docs: copulas#clayton-lower-tail-dependence
//! Backs the Clayton example on the copulas catalog page.

use stochastic_rs::copulas::bivariate::clayton::Clayton;
use stochastic_rs::copulas::correlation::kendall_tau;
use stochastic_rs::traits::BivariateExt;

#[test]
fn clayton_set_tau_then_sample() {
  let mut cop = Clayton::new();
  cop.set_tau(0.5); // tau => theta via Kendall inversion
  cop.set_theta(cop.compute_theta());
  let uv = cop.sample_with_seed(10_000, 42).unwrap(); // Array2<f64>, shape (10_000, 2)
  assert_eq!(uv.dim(), (10_000, 2));
  assert!(uv.iter().all(|&x| (0.0..=1.0).contains(&x)));

  let corr = kendall_tau(&uv); // pairwise Kendall tau matrix, shape (2, 2)
  let tau_hat = corr[[0, 1]];
  assert!((tau_hat - 0.5).abs() < 0.05, "tau_hat = {tau_hat}");
}

import stochastic_rs as srs

cop = srs.Clayton(theta=2.0)
uv = cop.sample(10_000, seed=42)              # shape (10_000, 2)
# kendall_tau_matrix takes an (n, k) matrix and returns the k×k pairwise-tau
# matrix; uv is already (n, 2), so [0, 1] is the one off-diagonal entry.
tau = srs.kendall_tau_matrix(uv)[0, 1]
print("Kendall's tau ≈", tau)                 # ≈ 0.5

// docs: copulas#bb1-bb7-two-parameter-tail-asymmetry
//! Backs the BB1 / vine-fitting example on the copulas catalog page.

use stochastic_rs::copulas::bivariate::bb1::Bb1;
use stochastic_rs::copulas::gof::gof_cramer_von_mises;
use stochastic_rs::copulas::gof::pseudo_observations;
use stochastic_rs::copulas::multivariate::fit::PairFamily;
use stochastic_rs::copulas::multivariate::fit::SelectionCriterion;
use stochastic_rs::copulas::multivariate::fit::VineStructure;
use stochastic_rs::copulas::multivariate::fit::fit_vine;
use stochastic_rs::traits::BivariateExt;

#[test]
fn bb1_fit_gof_and_vine_selection() {
  // BB1 with lower tail 2^{-1/(θδ)} and upper tail 2 - 2^{1/δ}: asymmetric tails.
  let truth = Bb1::new(Some(0.8), Some(1.6), None);
  let tails = truth.tail_dependence();
  assert!(tails.lower > 0.5 && tails.upper > 0.4 && tails.lower != tails.upper);
  // Every fit and test runs on pseudo-observations (normalised ranks), as with real data.
  let uv = pseudo_observations(&truth.sample_with_seed(2_000, 42).unwrap()); // shape (2_000, 2)

  // Maximum-likelihood fit of (θ, δ) recovers the generating pair.
  let mut fitted = Bb1::default();
  fitted.fit(&uv).unwrap();
  assert!((fitted.theta.unwrap() - 0.8).abs() < 0.25 && (fitted.delta - 1.6).abs() < 0.25);

  // Parametric-bootstrap Cramér–von Mises test does not reject the true family.
  let gof = gof_cramer_von_mises(&fitted, &uv, 20, 7, |c, x| c.fit(x)).unwrap();
  assert!(gof.p_value > 0.05, "p = {}", gof.p_value);

  // A two-column D-vine fit selects a family on the single edge by AIC.
  let fit = fit_vine(
    &uv,
    VineStructure::DVine,
    &PairFamily::ALL,
    SelectionCriterion::Aic,
  )
  .unwrap();
  assert_eq!(fit.families.len(), 1);
  assert!(fit.aic.is_finite() && fit.log_likelihood > 0.0);
}

import numpy as np
import stochastic_rs as srs

truth = srs.Bb1(theta=0.8, delta=1.6)
uv = srs.pseudo_observations(truth.sample(2_000, seed=42))
lower, upper = truth.tail_dependence()      # asymmetric: ≈ 0.58 vs ≈ 0.46

fitted = srs.Bb1()
fitted.fit(uv)                              # MLE of (theta, delta)
stat, p = srs.copula_gof("bb1", uv, replications=50, seed=7)
print(fitted.theta(), fitted.delta(), p)    # p well above 0.05

# 3-column data: structure + family selection with BIC
u3 = srs.pseudo_observations(np.column_stack([uv, np.random.default_rng(1).uniform(size=2_000)]))
fit = srs.fit_vine(u3, structure="dvine", criterion="bic")
print(fit["order"], fit["families"], fit["aic"])

// docs: copulas#gaussian-copula-multivariate-sampling
//! Backs the Gaussian multivariate copula example on the copulas catalog
//! page. The Cholesky factorisation runs on the pure-Rust `faer`, so the
//! file is gated on that feature.

use ndarray::array;
use stochastic_rs::copulas::multivariate::gaussian::GaussianMultivariate;
use stochastic_rs::traits::MultivariateExt;

#[test]
fn gaussian_multivariate_sample() {
  let corr = array![[1.0, 0.7, 0.3], [0.7, 1.0, 0.5], [0.3, 0.5, 1.0]];
  let cop = GaussianMultivariate::new_with_corr(corr).unwrap();
  let samples = cop.sample(10_000).unwrap(); // Array2<f64>, shape (10_000, 3)
  assert_eq!(samples.dim(), (10_000, 3));
}