### ----------------------------------------------------
### -- Math::Libgsl::RandomDistribution
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/01-raw-randist.rakutest
### ----------------------------------------------------
### -- Chunk 2 of 2
-tol(10⁻¹⁰), 'pareto with pdf';
  is-approx gsl_cdf_pareto_P(3e0, 1.9e0, 2.75e0), 0.15237894848417666e0, :abs-tol(10⁻¹⁰), 'P cumulative pareto';
  is-approx gsl_cdf_pareto_Q(3e0, 1.9e0, 2.75e0), 0.8476210515158233e0, :abs-tol(10⁻¹⁰), 'Q cumulative pareto';
  is-approx gsl_cdf_pareto_Pinv(0.15237894848417666e0, 1.9e0, 2.75e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative pareto';
  is-approx gsl_cdf_pareto_Qinv(0.8476210515158233e0, 1.9e0, 2.75e0), 3e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative pareto';
  gsl_rng_free($r);
}

subtest 'spherical vector distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  my num64 ($x, $y) = (0e0, 0e0);
  gsl_ran_dir_2d($r, $x, $y);
  ok $x ≅ -0.617745613497854e0 && $y ≅ -0.7863779988047479e0, '2D random vector';

  ($x, $y) = (0e0, 0e0);
  gsl_ran_dir_2d_trig_method($r, $x, $y);
  ok $x ≅ 0.11500033916619544e0 && $y ≅ 0.9933654523848008e0, '2D random vector trig method';

  my num64 $z;
  ($x, $y, $z) = (0e0, 0e0, 0e0);
  gsl_ran_dir_3d($r, $x, $y, $z);
  ok $x ≅ -0.024188741233946875e0 && $y ≅ 0.7364240456943486e0 && $z ≅ 0.6760876642275653e0, '3D random vector';

  my $n = 5;
  my $xarr = CArray[num64].new: 0e0 xx $n;
  gsl_ran_dir_nd($r, $n, $xarr);
  ok ([&&] ($xarr[^$n].list
    Z≅
    (0.1717392856031707e0, 0.547345960645299e0, -0.814904936760339e0, 0.05741225862147077e0, -0.059596927348619114e0))),
    'nD random vector';

  gsl_rng_free($r);
}

subtest 'Weibull distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_weibull($r, 3.14e0, 2.75e0), 0.15568189304354582e0, :abs-tol(10⁻¹⁰), 'weibull';
  is-approx gsl_ran_weibull_pdf(2e0, 3.14e0, 2.75e0), 0.2978229948509798e0, :abs-tol(10⁻¹⁰), 'weibull with pdf';
  is-approx gsl_cdf_weibull_P(2e0, 3.14e0, 2.75e0), 0.25117631509661126e0, :abs-tol(10⁻¹⁰), 'P cumulative weibull';
  is-approx gsl_cdf_weibull_Q(2e0, 3.14e0, 2.75e0), 0.7488236849033888e0, :abs-tol(10⁻¹⁰), 'Q cumulative weibull';
  is-approx gsl_cdf_weibull_Pinv(0.25117631509661126e0, 3.14e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative weibull';
  is-approx gsl_cdf_weibull_Qinv(0.7488236849033888e0, 3.14e0, 2.75e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative weibull';
  gsl_rng_free($r);
}

subtest 'type-1 Gumbel distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gumbel1($r, 3.12e0, 4.56e0), 3.134221698863258e0, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gsl_ran_gumbel1_pdf(2e0, 3.12e0, 4.56e0), 0.02749542323890301e0, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gsl_cdf_gumbel1_P(2e0, 3.12e0, 4.56e0), 0.9911480698975662e0, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Q(2e0, 3.12e0, 4.56e0), 0.008851930102433793e0, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Pinv(0.9911480698975662e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gsl_cdf_gumbel1_Qinv(0.008851930102433793e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel1';
  gsl_rng_free($r);
}

subtest 'type-2 Gumbel distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  is-approx gsl_ran_gumbel2($r, 3.12e0, 4.56e0), 22.970750721534436e0, :abs-tol(10⁻¹⁰), 'gumbel1';
  is-approx gsl_ran_gumbel2_pdf(2e0, 3.12e0, 4.56e0), 0.4842675579059279e0, :abs-tol(10⁻¹⁰), 'gumbel1 with pdf';
  is-approx gsl_cdf_gumbel2_P(2e0, 3.12e0, 4.56e0), 0.5918470962288894e0, :abs-tol(10⁻¹⁰), 'P cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Q(2e0, 3.12e0, 4.56e0), 0.4081529037711105e0, :abs-tol(10⁻¹⁰), 'Q cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Pinv(0.5918470962288894e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse P cumulative gumbel1';
  is-approx gsl_cdf_gumbel2_Qinv(0.4081529037711105e0, 3.12e0, 4.56e0), 2e0, :abs-tol(10⁻¹⁰), 'inverse Q cumulative gumbel2';
  gsl_rng_free($r);
}

subtest 'Dirichlet distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $K = 2;
  my $alpha = CArray[num64].new: 2.5e0, 5e0;
  my $theta = CArray[num64].new: 0e0, 0e0;

  gsl_ran_dirichlet($r, $K, $alpha, $theta);
  ok ([&&] ($theta[^$K].list Z≅ (0.37979184993320686e0, 0.6202081500667932e0))), 'dirichlet';

  $theta[0] = .2e0;
  $theta[1] = .8e0;
  ok gsl_ran_dirichlet_pdf($K, $alpha, $theta) ≅ 2.148771883658201e0, 'dirichlet with pdf';

  ok gsl_ran_dirichlet_lnpdf($K, $alpha, $theta) ≅ 0.7648964620298799e0, 'ln dirichlet with pdf';

  gsl_rng_free($r);
}

subtest 'general discrete distribution' => {
  my $*TOLERANCE = 10⁻¹⁰;
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $P = CArray[num64].new: .59e0, .4e0, .01e0;

  my $g = gsl_ran_discrete_preproc(3, $P);
  isa-ok $g, gsl_ran_discrete_t, 'create lookup table';
  ok gsl_ran_discrete($r, $g) ~~ 0|1|2, 'get discrete random number';
  ok gsl_ran_discrete_pdf(2, $g) ≅ .01e0, 'probability of a value';
  lives-ok { gsl_ran_discrete_free($g) }, 'destroy lookup table ';

  gsl_rng_free($r);
}

subtest 'Poisson distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_poisson($r, 30e0) == 34, 'poisson';
  is-approx gsl_ran_poisson_pdf(2, 5e0), 0.08422433748856832e0, :abs-tol(10⁻¹⁰), 'poisson with pdf';
  is-approx gsl_cdf_poisson_P(2, 5e0), 0.12465201948308077e0, :abs-tol(10⁻¹⁰), 'P cumulative poisson';
  is-approx gsl_cdf_poisson_Q(2, 5e0), 0.8753479805169192e0, :abs-tol(10⁻¹⁰), 'Q cumulative poisson';

  gsl_rng_free($r);
}

subtest 'Bernoulli distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_bernoulli($r, .3e0) == 0, 'bernoulli';
  is-approx gsl_ran_bernoulli_pdf(1, .3e0), 0.3e0, :abs-tol(10⁻¹⁰), 'bernoulli with pdf';

  gsl_rng_free($r);
}

subtest 'binomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_binomial($r, .3e0, 5) == 5, 'binomial';
  is-approx gsl_ran_binomial_pdf(2, .3e0, 5), 0.3086999999999999e0, :abs-tol(10⁻¹⁰), 'binomial with pdf';
  is-approx gsl_cdf_binomial_P(2, .3e0, 5), 0.83692e0, :abs-tol(10⁻¹⁰), 'P cumulative binomial';
  is-approx gsl_cdf_binomial_Q(2, .3e0, 5), 0.16307999999999997e0, :abs-tol(10⁻¹⁰), 'Q cumulative binomial';

  gsl_rng_free($r);
}

subtest 'multinomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my $K = 3;
  my $sum_n = 100;
  my $p = CArray[num64].new: 2e0, 7e0, 1e0;
  my $n = CArray[uint32].new: 0 xx 3;

  gsl_ran_multinomial($r, $K, $sum_n, $p, $n);
  is-deeply $n.list, (20, 72, 8), 'multinomial';
  is-approx gsl_ran_multinomial_pdf($K, $p, $n), 0.011455694618856253, :abs-tol(10⁻¹⁰), 'multinomial with pdf';
  is-approx gsl_ran_multinomial_lnpdf($K, $p, $n), -4.469268325992729, :abs-tol(10⁻¹⁰), 'ln multinomial with pdf';

  gsl_rng_free($r);
}

subtest 'negative binomial distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_negative_binomial($r, .3e0, 20e0) == 45, 'negative binomial';
  is-approx gsl_ran_negative_binomial_pdf(20, .3e0, 20e0), 0.0019175722376507298e0, :abs-tol(10⁻¹⁰), 'negative binomial with pdf';
  is-approx gsl_cdf_negative_binomial_P(20, .3e0, 20e0), 0.006254504372434949e0, :abs-tol(10⁻¹⁰), 'P cumulative negative binomial';
  is-approx gsl_cdf_negative_binomial_Q(20, .3e0, 20e0), 0.9937454956275651e0, :abs-tol(10⁻¹⁰), 'Q cumulative negative binomial';

  gsl_rng_free($r);
}

subtest 'Pascal distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_pascal($r, .8e0, 3) == 0, 'pascal';
  is-approx gsl_ran_pascal_pdf(2, .8e0, 3), 0.12287999999999971e0, :abs-tol(10⁻¹⁰), 'pascal with pdf';
  is-approx gsl_cdf_pascal_P(2, .8e0, 3), 0.94208e0, :abs-tol(10⁻¹⁰), 'P cumulative pascal';
  is-approx gsl_cdf_pascal_Q(2, .8e0, 3), 0.057920000000000006e0, :abs-tol(10⁻¹⁰), 'Q cumulative pascal';

  gsl_rng_free($r);
}

subtest 'geometric distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_geometric($r, .5e0) == 1, 'geometric';
  is-approx gsl_ran_geometric_pdf(2, .5e0), .25e0, :abs-tol(10⁻¹⁰), 'geometric with pdf';
  is-approx gsl_cdf_geometric_P(2, .5e0), .75e0, :abs-tol(10⁻¹⁰), 'P cumulative geometric';
  is-approx gsl_cdf_geometric_Q(2, .5e0), .25e0, :abs-tol(10⁻¹⁰), 'Q cumulative geometric';

  gsl_rng_free($r);
}

subtest 'hypergeometric distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_hypergeometric($r, 5, 7, 4) == 2, 'hypergeometric';
  is-approx gsl_ran_hypergeometric_pdf(2, 5, 7, 4), 0.4242424242424253e0, :abs-tol(10⁻¹⁰), 'hypergeometric with pdf';
  is-approx gsl_cdf_hypergeometric_P(2, 5, 7, 4), 0.8484848484848483e0, :abs-tol(10⁻¹⁰), 'P cumulative hypergeometric';
  is-approx gsl_cdf_hypergeometric_Q(2, 5, 7, 4), 0.1515151515151517e0, :abs-tol(10⁻¹⁰), 'Q cumulative hypergeometric';

  gsl_rng_free($r);
}

subtest 'logarithmic distribution' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);

  ok gsl_ran_logarithmic($r, .4e0) == 1, 'logarithmic';
  is-approx gsl_ran_logarithmic_pdf(2, .4e0), 0.15660921511769743e0, :abs-tol(10⁻¹⁰), 'logarithmic with pdf';

  gsl_rng_free($r);
}

subtest 'Wishart distribution' => {
  my gsl_matrix $V    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $L    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $X    = gsl_matrix_calloc(2, 2);
  my gsl_matrix $L_X  = gsl_matrix_calloc(2, 2);
  my gsl_matrix $work = gsl_matrix_calloc(2, 2);
  my num64 $obs_res;
  my num64 $df = 3e0;
  my $*TOLERANCE = 10⁻¹⁰;

  gsl_matrix_set($V, 0, 0, 1e0);
  gsl_matrix_set($V, 1, 1, 1e0);
  gsl_matrix_set($V, 0, 1, .3e0);
  gsl_matrix_set($V, 1, 0, .3e0);

  gsl_matrix_memcpy($L, $V);
  gsl_linalg_cholesky_decomp1($L);

  gsl_matrix_set($X, 0, 0, 2.213322e0);
  gsl_matrix_set($X, 1, 1, 3.285779e0);
  gsl_matrix_set($X, 0, 1, 1.453357e0);
  gsl_matrix_set($X, 1, 0, 1.453357e0);

  gsl_matrix_memcpy($L_X, $X);
  gsl_linalg_cholesky_decomp1($L_X);

  gsl_ran_wishart_log_pdf($X, $L_X, $df, $L, $obs_res, $work);
  ok $obs_res ≅ -4.931913612377813e0, 'log wishart with pdf';

  gsl_ran_wishart_pdf($X, $L_X, $df, $L, $obs_res, $work);
  ok $obs_res ≅ 0.007212687778224e0, 'wishart with pdf';

  gsl_matrix_free($V);
  gsl_matrix_free($L);
  gsl_matrix_free($X);
  gsl_matrix_free($L_X);
  gsl_matrix_free($work);
}

subtest 'shuffling and sampling' => {
  my gsl_rng $r = mgsl_rng_setup(DEFAULT);
  my CArray[int32] $x .= new: ^10;
  my CArray[int32] $y .= new: ^3;

  gsl_ran_shuffle($r, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $x[^10], (5, 7, 2, 3, 4, 1, 0, 6, 9, 8), 'shuffle';

  gsl_ran_choose($r, pointer-to($y), 3, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $y[^3], (2, 3, 9), 'choose';

  gsl_ran_sample($r, pointer-to($y), 3, pointer-to($x), 10, nativesizeof(int32).Int) for ^10;
  is-deeply $y[^3], (3, 3, 8), 'sample';

  gsl_rng_free($r);
}

done-testing;


