### ----------------------------------------------------
### -- Math::Libgsl::RandomDistribution
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-randist.rakutest
### ----------------------------------------------------
### -- Chunk 1 of 2
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::Random;
use Math::Libgsl::Matrix;
use Math::Libgsl::BLAS;
use Math::Libgsl::LinearAlgebra;
use Math::Libgsl::Constants;
use Math::Libgsl::RandomDistribution;

subtest 'gaussian distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gaussian($r, 3), 0.4017558243560277, :rel-tol(10⁻¹²), 'gaussian random variate';
  is-approx gaussian-pdf(1, 3), 0.12579440923099774, :rel-tol(10⁻¹²), 'gaussian probability density';
  is-approx gaussian-ziggurat($r, 3), 1.3570084962645521, :rel-tol(10⁻¹²), 'gaussian ziggurat';
  is-approx gaussian-ratio-method($r, 3), 4.571692762640136, :rel-tol(10⁻¹²), 'gaussian ratio method';
  is-approx ugaussian($r), 0.2712258328755056, :rel-tol(10⁻¹²), 'unit gaussian';
  is-approx ugaussian-pdf(3), 0.0044318484119380075, :rel-tol(10⁻¹²), 'unit gaussian probability density';
  is-approx ugaussian-ratio-method($r), 1.7149190276802055, :rel-tol(10⁻¹²), 'unit gaussian ratio method';
  is-approx gaussian-P(1, 3), 0.6305586598182363, :rel-tol(10⁻¹²), 'cumulative distribution P(x)';
  is-approx gaussian-Q(1, 3), 0.36944134018176367, :rel-tol(10⁻¹²), 'cumulative distribution Q(x)';
  is-approx gaussian-Pinv(0.6305586598182363, 3), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution P(x)';
  is-approx gaussian-Qinv(0.36944134018176367, 3), 1, :rel-tol(10⁻¹²), 'inverse cumulative distribution Q(x)';
  is-approx ugaussian-P(1), 0.8413447460685429, :rel-tol(10⁻¹²), 'unit cumulative distribution P(x)';
  is-approx ugaussian-Q(1), 0.15865525393145705, :rel-tol(10⁻¹²), 'unit cumulative distribution Q(x)';
  is-approx ugaussian-Pinv(0.8413447460685429), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution P(x)';
  is-approx ugaussian-Qinv(0.15865525393145705), 1, :rel-tol(10⁻¹²), 'inverse unit cumulative distribution Q(x)';
}

subtest 'gaussian tail distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gaussian-tail($r, 1, 3), 5.023225218761322, :rel-tol(10⁻¹²), 'gaussian tail';
  is-approx gaussian-tail-pdf(1, 1, 3), 0.3404990063350988, :rel-tol(10⁻¹²), 'gaussian tail probability density';
  is-approx ugaussian-tail($r, 1), 1.198044133837379, :rel-tol(10⁻¹²), 'unit gaussian tail';
  is-approx ugaussian-tail-pdf(1, 1), 1.525135276160981, :rel-tol(10⁻¹²), 'unit gaussian tail probability density';
}

subtest 'bivariate gaussian distribution' => {
  my $*TOLERANCE = 10⁻¹²;
  my Math::Libgsl::Random $r .= new;
  my ($x, $y) = bivariate-gaussian($r, 3, 2, .3);
  ok ([&&] (($x, $y) Z≅ (-0.1952914837346669, 0.21644212342026375))), 'bivariate gaussian';
  is-approx bivariate-gaussian-pdf(0e0, 0e0, 3e0, 2e0, .3e0), 0.027806618922095617, :rel-tol(10⁻¹²), 'bivariate gaussian probability density';
}

subtest 'bivariate gaussian distribution' => {
  my $*TOLERANCE = 10⁻¹²;
  my $N = 100;
  my $d = 2;
  my Math::Libgsl::Random $r .= new;
  my Math::Libgsl::Matrix $samples .= new: $N, $d;
  my Math::Libgsl::Matrix $L       .= new: $d, $d;
  my Math::Libgsl::Vector $mu      .= new: $d;
  my Math::Libgsl::Vector $x       .= new: $d;
  $mu[0] = 1; $mu[1] = 2;
  $L[0;0] = 4; $L[0;1] = 2; $L[1;0] = 2; $L[1;1] = 3;
  cholesky-decomp1($L);

  for ^$N -> $i {
    my $sample = multivariate-gaussian($r, $mu, $L);
    $samples.set-row($i, $sample);
  }
  ok ([&&] ($samples.get-row(0) Z≅ (1.2678372162373517, 2.0093249906121344))), 'multivariate gaussian';
  my $mu-hat = multivariate-gaussian-mean($samples);
  my $L-hat  = multivariate-gaussian-vcov($samples);
  ok ([&&] ($mu-hat[^2] Z≅ (1.2705286868498327, 2.2420437160261))), 'multivariate gaussian mean';
  ok ([&&] ((gather for ^2 X ^2 -> ($i, $j) { take $L-hat[$i;$j] })
      Z≅
      (4.105455552558069, 2.097351009513896, 2.097351009513896, 3.229457820732376))), 'multivariate gaussian vcov';
  is-approx multivariate-gaussian-pdf($x, $mu, $L), 0.028294217120391e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';
  is-approx multivariate-gaussian-log-pdf($x, $mu, $L), -3.565097837249263e0, :abs-tol(10⁻¹⁰), 'multivariate gaussian with pdf';
}

subtest 'exponential distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx exponential($r, 2), 16.523156432740787, :abs-tol(10⁻¹⁰), 'exponential';
  is-approx exponential-pdf(2, 2), 0.18393972058572117, :abs-tol(10⁻¹⁰), 'exponential with pdf';
  is-approx exponential-P(2, 2), 0.6321205588285577, :abs-tol(10⁻¹⁰), 'P cumulative exponential';
  is-approx exponential-Q(2, 2), 0.36787944117144233, :abs-tol(10⁻¹⁰), 'Q cumulative exponential';
  is-approx exponential-Pinv(0.6321205588285577, 2), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative exponential';
  is-approx exponential-Qinv(0.36787944117144233, 2), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative exponential';
}

subtest 'Laplace distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx laplace($r, 2.75), 0.001420747954609721, :abs-tol(10⁻¹⁰), 'laplace';
  is-approx laplace-pdf(2, 2.75), 0.08785910567087736, :abs-tol(10⁻¹⁰), 'laplace with pdf';
  is-approx laplace-P(2, 2.75), 0.7583874594050872, :abs-tol(10⁻¹⁰), 'P cumulative laplace';
  is-approx laplace-Q(2, 2.75), 0.24161254059491272, :abs-tol(10⁻¹⁰), 'Q cumulative laplace';
  is-approx laplace-Pinv(0.7583874594050872, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative laplace';
  is-approx laplace-Qinv(0.24161254059491272, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative laplace';
}

subtest 'exponential power distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx exppow($r, 3.7, 0.3), 96.18418180115754, :abs-tol(10⁻¹⁰), 'exppow';
  is-approx exppow-pdf(2, 3.7, 0.3), 0.006353720526045815, :abs-tol(10⁻¹⁰), 'exppow with pdf';
  is-approx exppow-P(2, 3.7, 0.3), 0.5155820414218526, :abs-tol(10⁻¹⁰), 'P cumulative expow';
  is-approx exppow-Q(2, 3.7, 0.3), 0.4844179585781474, :abs-tol(10⁻¹⁰), 'Q cumulative expow';
}

subtest 'Cauchy distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx cauchy($r, 2), -0.001622639831190868, :abs-tol(10⁻¹⁰), 'cauchy';
  is-approx cauchy-pdf(2, 2), 0.07957747154594767, :abs-tol(10⁻¹⁰), 'cauchy with pdf';
  is-approx cauchy-P(2, 2), 0.75, :abs-tol(10⁻¹⁰), 'P cumulative cauchy';
  is-approx cauchy-Q(2, 2), 0.25, :abs-tol(10⁻¹⁰), 'Q cumulative cauchy';
  is-approx cauchy-Pinv(.75, 2), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative cauchy';
  is-approx cauchy-Qinv(.25, 2), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative cauchy';
}

subtest 'Rayleigh distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx rayleigh($r, 1.9), 0.04318348873449323, :abs-tol(10⁻¹⁰), 'rayleigh';
  is-approx rayleigh-pdf(2, 1.9), 0.3183584873991547, :abs-tol(10⁻¹⁰), 'rayleigh with pdf';
  is-approx rayleigh-P(2, 1.9), 0.4253629302445258, :abs-tol(10⁻¹⁰), 'P cumulative rayleigh';
  is-approx rayleigh-Q(2, 1.9), 0.5746370697554741, :abs-tol(10⁻¹⁰), 'Q cumulative rayleigh';
  is-approx rayleigh-Pinv(0.4253629302445258, 1.9), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative rayleigh';
  is-approx rayleigh-Qinv(0.5746370697554741, 1.9), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative rayleigh';
}

subtest 'Rayleigh tail distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx rayleigh-tail($r, 2.7, 1.9), 2.7003453137884574, :abs-tol(10⁻¹⁰), 'rayleigh tail';
  is-approx rayleigh-tail-pdf(2, 2.7, 1.9), 0, :abs-tol(10⁻¹⁰), 'rayleigh tail with pdf';
}

subtest 'Landau distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx landau($r), 3880.037426254597, :abs-tol(10⁻¹⁰), 'landau';
  is-approx landau-pdf(2), 0.10491298949080555, :abs-tol(10⁻¹⁰), 'landau with pdf';
}

subtest 'Levy alpha-Stable distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx levy($r, 5, 1), 6162.79707164862, :abs-tol(10⁻¹⁰), 'levy';
}

subtest 'Levy skew alpha-Stable distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx levy-skew($r, 5, 1, 0), 6162.79707164862, :abs-tol(10⁻¹⁰), 'levy skew';
}

subtest 'gamma distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gamma($r, 2.5, 2.17), 7.592395030767021, :abs-tol(10⁻¹⁰), 'gamma';
  is-approx gamma-knuth($r, 2.5, 2.17), 7.554777177884447, :abs-tol(10⁻¹⁰), 'gamma knuth';
  is-approx gamma-pdf(2, 2.5, 2.17), 0.12203602332361416, :abs-tol(10⁻¹⁰), 'gamma with pdf';
  is-approx gamma-P(2, 2.5, 2.17), 0.12962769871792748, :abs-tol(10⁻¹⁰), 'P cumulative gamma';
  is-approx gamma-Q(2, 2.5, 2.17), 0.8703723012820725, :abs-tol(10⁻¹⁰), 'Q cumulative gamma';
  is-approx gamma-Pinv(0.12962769871792748, 2.5, 2.17), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative gamma';
  is-approx gamma-Qinv(0.8703723012820725, 2.5, 2.17), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gamma';
}

subtest 'flat distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx flat($r, 3, 4), 3.999741748906672, :abs-tol(10⁻¹⁰), 'flat';
  is-approx flat-pdf(2, 3, 4), 0, :abs-tol(10⁻¹⁰), 'flat with pdf';
  is-approx flat-P(2, 3, 4), 0, :abs-tol(10⁻¹⁰), 'P cumulative flat';
  is-approx flat-Q(2, 3, 4), 1, :abs-tol(10⁻¹⁰), 'Q cumulative flat';
  is-approx flat-Pinv(0, 3, 4), 3, :abs-tol(10⁻¹⁰), 'inverse P cumulative flat';
  is-approx flat-Qinv(1, 3, 4), 3, :abs-tol(10⁻¹⁰), 'inverse Q cumulative flat';
}

subtest 'lognormal distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx lognormal($r, 2.7, 3), 12.239990687580857, :abs-tol(10⁻¹⁰), 'lognormal';
  is-approx lognormal-pdf(2, 2.7, 3), 0.053160178764929525, :abs-tol(10⁻¹⁰), 'lognormal with pdf';
  is-approx lognormal-P(2, 2.7, 3), 0.25176338701440004, :abs-tol(10⁻¹⁰), 'P cumulative lognormal';
  is-approx lognormal-Q(2, 2.7, 3), 0.7482366129856, :abs-tol(10⁻¹⁰), 'Q cumulative lognormal';
  is-approx lognormal-Pinv(0.25176338701440004, 2.7, 3), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative lognormal';
  is-approx lognormal-Qinv(0.7482366129856, 2.7, 3), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative lognormal';
}

subtest 'chi-squared distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx chisq($r, 13), 16.53548783112717, :abs-tol(10⁻¹⁰), 'chisq';
  is-approx chisq-pdf(2, 13), 0.0006389340989638759, :abs-tol(10⁻¹⁰), 'chisq with pdf';
  is-approx chisq-P(2, 13), 0.00022625008465604695, :abs-tol(10⁻¹⁰), 'P cumulative chisq';
  is-approx chisq-Q(2, 13), 0.999773749915344, :abs-tol(10⁻¹⁰), 'Q cumulative chisq';
  is-approx chisq-Pinv(0.00022625008465604695, 13), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative chisq';
  is-approx chisq-Qinv(0.999773749915344, 13), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative chisq';
}

subtest 'f-distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx fdist($r, 3, 4), 1.265383001140288, :abs-tol(10⁻¹⁰), 'fdist';
  is-approx fdist-pdf(2, 3, 4), 0.13942740046346694, :abs-tol(10⁻¹⁰), 'fdist with pdf';
  is-approx fdist-P(2, 3, 4), 0.743612802471824, :abs-tol(10⁻¹⁰), 'P cumulative fdist';
  is-approx fdist-Q(2, 3, 4), 0.2563871975281759, :abs-tol(10⁻¹⁰), 'Q cumulative fdist';
  is-approx fdist-Pinv(0.743612802471824, 3, 4), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative fdist';
  is-approx fdist-Qinv(0.2563871975281759, 3, 4), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative fdist';
}

subtest 't-distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx tdist($r, 1.75), 0.2825086397013947, :abs-tol(10⁻¹⁰), 'tdist';
  is-approx tdist-pdf(2, 1.75), 0.06778157292044809, :abs-tol(10⁻¹⁰), 'tdist with pdf';
  is-approx tdist-P(2, 1.75), 0.8993101012459747, :abs-tol(10⁻¹⁰), 'P cumulative tdist';
  is-approx tdist-Q(2, 1.75), 0.1006898987540253, :abs-tol(10⁻¹⁰), 'Q cumulative tdist';
  is-approx tdist-Pinv(0.8993101012459747, 1.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative tdist';
  is-approx tdist-Qinv(0.1006898987540253, 1.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative tdist';
}

subtest 'beta distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx beta($r, 2, 3), 0.45158505658193976, :abs-tol(10⁻¹⁰), 'beta';
  is-approx beta-pdf(2, 2, 3), 0, :abs-tol(10⁻¹⁰), 'beta with pdf';
  is-approx beta-P(2, 2, 3), 1, :abs-tol(10⁻¹⁰), 'P cumulative beta';
  is-approx beta-Q(2, 2, 3), 0, :abs-tol(10⁻¹⁰), 'Q cumulative beta';
  is-approx beta-Pinv(1, 2, 3), 1, :abs-tol(10⁻¹⁰), 'inverse P cumulative beta';
  is-approx beta-Qinv(0, 2, 3), 1, :abs-tol(10⁻¹⁰), 'inverse Q cumulative beta';
}

subtest 'logistic distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx logistic($r, 3.1), 25.61009178896598, :abs-tol(10⁻¹⁰), 'logistic';
  is-approx logistic-pdf(2, 3.1), 0.07280296434630028, :abs-tol(10⁻¹⁰), 'logistic with pdf';
  is-approx logistic-P(2, 3.1), 0.6559192436053651, :abs-tol(10⁻¹⁰), 'P cumulative logistic';
  is-approx logistic-Q(2, 3.1), 0.344080756394635, :abs-tol(10⁻¹⁰), 'Q cumulative logistic';
  is-approx logistic-Pinv(0.6559192436053651, 3.1), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative logistic';
  is-approx logistic-Qinv(0.344080756394635, 3.1), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative logistic';
}

subtest 'Pareto distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx pareto($r, 1.9, 2.75), 2.7503738581610317, :abs-tol(10⁻¹⁰), 'pareto';
  is-approx pareto-pdf(3, 1.9, 2.75), 0.5368266659600216, :abs-tol(10⁻¹⁰), 'pareto with pdf';
  is-approx pareto-P(3, 1.9, 2.75), 0.15237894848417666, :abs-tol(10⁻¹⁰), 'P cumulative pareto';
  is-approx pareto-Q(3, 1.9, 2.75), 0.8476210515158233, :abs-tol(10⁻¹⁰), 'Q cumulative pareto';
  is-approx pareto-Pinv(0.15237894848417666, 1.9, 2.75), 3, :abs-tol(10⁻¹⁰), 'inverse P cumulative pareto';
  is-approx pareto-Qinv(0.8476210515158233, 1.9, 2.75), 3, :abs-tol(10⁻¹⁰), 'inverse Q cumulative pareto';
}

subtest 'spherical vector distribution' => {
  my Math::Libgsl::Random $r .= new;
  my $*TOLERANCE = 10⁻¹⁰;
  my ($x, $y) = d2dir($r);
  ok $x ≅ -0.617745613497854 && $y ≅ -0.7863779988047479, '2D random vector';
  ($x, $y) = d2dir-trig-method($r);
  ok $x ≅ 0.11500033916619544 && $y ≅ 0.9933654523848008, '2D random vector trig method';
  my $z;
  ($x, $y, $z) = d3dir($r);
  ok $x ≅ -0.024188741233946875 && $y ≅ 0.7364240456943486 && $z ≅ 0.6760876642275653, '3D random vector';
  my @ret = dndir($r, 5);
  ok ([&&] (@ret[^5]
    Z≅
    (0.1717392856031707, 0.547345960645299, -0.814904936760339, 0.05741225862147077, -0.059596927348619114))),
    'nD random vector';
}

subtest 'Weibull distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx weibull($r, 3.14, 2.75), 0.15568189304354582, :abs-tol(10⁻¹⁰), 'weibull';
  is-approx weibull-pdf(2, 3.14, 2.75), 0.2978229948509798, :abs-tol(10⁻¹⁰), 'weibull with pdf';
  is-approx weibull-P(2, 3.14, 2.75), 0.25117631509661126, :abs-tol(10⁻¹⁰), 'P cumulative weibull';
  is-approx weibull-Q(2, 3.14, 2.75), 0.7488236849033888, :abs-tol(10⁻¹⁰), 'Q cumulative weibull';
  is-approx weibull-Pinv(0.25117631509661126, 3.14, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative weibull';
  is-approx weibull-Qinv(0.7488236849033888, 3.14, 2.75), 2, :abs-tol(10⁻¹⁰), 'inverse Q cumulative weibull';
}

subtest 'type-1 Gumbel distribution' => {
  my Math::Libgsl::Random $r .= new;
  is-approx gumbel1($r, 3.12, 4.56), 3.134221698863258, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gumbel1-pdf(2, 3.12, 4.56), 0.02749542323890301, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gumbel1-P(2, 3.12, 4.56), 0.9911480698975662, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gumbel1-Q(2, 3.12, 4.56), 0.008851930102433793, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gumbel1-Pinv(0.9911480698975662, 3.12, 4.56), 2, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gumbel1