### ----------------------------------------------------
### -- Math::Libgsl::Wavelet
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/02-raw-2d.rakutest
### ----------------------------------------------------
#!/usr/bin/env raku

use Test;
use lib 'lib';
use NativeCall;
use Math::Libgsl::Raw::Wavelet;
use Math::Libgsl::Raw::Matrix :ALL;
use Math::Libgsl::Constants;

sub test-num64-matrix(gsl_matrix $m1, gsl_matrix $m2 where $m1.size1 == $m2.size1 && $m1.size2 == $m2.size2 --> Bool)
{
  my Bool $result = True;
  my $*TOLERANCE = 10⁻¹⁰;
  for ^$m1.size1 X ^$m1.size2 -> ($i, $j) {
    $result &&= gsl_matrix_get($m1, $i, $j) ≅ gsl_matrix_get($m2, $i, $j);
  }
  $result
}

constant \N    = 256;
constant \kind =   4;
my $*TOLERANCE = 10⁻¹⁰;

my gsl_wavelet $w = mgsl_wavelet_setup(DAUBECHIES, kind);
my gsl_wavelet_workspace $work = gsl_wavelet_workspace_alloc(N);

subtest 'standard form on array' => {
  my $odata  = CArray[num64].allocate(N * N);
  my $idata = CArray[num64].allocate(N * N);

  for ^N X ^N -> ($i, $j) {
    $odata[$i * N + $j] = $idata[$i * N + $j] = ($i * N + $j).Num / (N * N);
  }

  ok gsl_wavelet2d_transform_forward($w, $odata, N, N, N, $work) == GSL_SUCCESS, 'forward transform';
  ok gsl_wavelet2d_transform_inverse($w, $odata, N, N, N, $work) == GSL_SUCCESS, 'inverse transform';
  ok ([&&] ($odata[^(N * N)] Z≅ $idata[^(N * N)])), 'comparison with original values';
}

subtest 'standard form on matrix' => {
  my gsl_matrix $idata = gsl_matrix_alloc(N, N);
  my gsl_matrix $odata = gsl_matrix_alloc(N, N);

  for ^N X ^N -> ($i, $j) {
    gsl_matrix_set($idata, $i, $j, ($i * N + $j).Num / (N * N));
    gsl_matrix_set($odata, $i, $j, ($i * N + $j).Num / (N * N));
  }

  ok gsl_wavelet2d_transform_matrix_forward($w, $odata, $work) == GSL_SUCCESS, 'forward transform';
  ok gsl_wavelet2d_transform_matrix_inverse($w, $odata, $work) == GSL_SUCCESS, 'inverse transform';
  ok test-num64-matrix($idata, $odata), 'comparison with original values';

  gsl_matrix_free($idata);
  gsl_matrix_free($odata);
}

subtest 'non standard form on array' => {
  my $odata  = CArray[num64].allocate(N * N);
  my $idata = CArray[num64].allocate(N * N);

  for ^N X ^N -> ($i, $j) {
    $odata[$i * N + $j] = $idata[$i * N + $j] = ($i * N + $j).Num / (N * N);
  }

  ok gsl_wavelet2d_nstransform_forward($w, $odata, N, N, N, $work) == GSL_SUCCESS, 'forward transform';
  ok gsl_wavelet2d_nstransform_inverse($w, $odata, N, N, N, $work) == GSL_SUCCESS, 'inverse transform';
  ok ([&&] ($odata[^(N * N)] Z≅ $idata[^(N * N)])), 'comparison with original values';
}

subtest 'non standard form on matrix' => {
  my gsl_matrix $idata = gsl_matrix_alloc(N, N);
  my gsl_matrix $odata = gsl_matrix_alloc(N, N);

  for ^N X ^N -> ($i, $j) {
    gsl_matrix_set($idata, $i, $j, ($i * N + $j).Num / (N * N));
    gsl_matrix_set($odata, $i, $j, ($i * N + $j).Num / (N * N));
  }

  ok gsl_wavelet2d_nstransform_matrix_forward($w, $odata, $work) == GSL_SUCCESS, 'forward transform';
  ok gsl_wavelet2d_nstransform_matrix_inverse($w, $odata, $work) == GSL_SUCCESS, 'inverse transform';
  ok test-num64-matrix($idata, $odata), 'comparison with original values';

  gsl_matrix_free($idata);
  gsl_matrix_free($odata);
}

gsl_wavelet_free($w);
gsl_wavelet_workspace_free($work);

done-testing;


