### ----------------------------------------------------
### -- Math::Libgsl::Matrix
### -- Licenses: Artistic-2.0
### -- Authors: Fernando Santagata
### -- File: t/04-matrix.rakutest
### ----------------------------------------------------
#!/usr/bin/env raku

use Test;
use lib 'lib';
use Math::Libgsl::Matrix;
use Math::Libgsl::Vector;
use Math::Libgsl::Constants;

subtest 'matrix allocation' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  isa-ok $m1, Math::Libgsl::Matrix, 'new using pairs';
  ok $m1.matrix.size1 == 3 && $m1.matrix.size2 == 3, 'size directly';
  is-deeply $m1.size, (3, 3), 'both sizes through method';
  ok $m1.size1 == 3, 'size1';
  ok $m1.size2 == 3, 'size2';
  my Math::Libgsl::Matrix $m2 .= new(3, 3);
  isa-ok $m2, Math::Libgsl::Matrix, 'new using values';
};

subtest 'matrix set & get' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  ok $m1.get(0, 0) == 0, 'get element';
  dies-ok { $m1.get(3, 3) }, 'get element (4,4) of a (3,3) matrix dies';
  ok $m1[0;0] == 0, 'get element using postcircumfix';
  dies-ok { $m1[3;3] }, 'get element (4,4) using postcircumfix of a (3,3) matrix dies';
  lives-ok { $m1.set(0, 0, 10) }, 'element set lives';
  ok $m1.get(0, 0) == 10e0, 'element set';
  dies-ok { $m1.set(3, 3, 10) }, 'set element (4,4) of a (3,3) matrix dies';
  lives-ok { $m1[0;1] = 12 }, 'set element using postcircumfix lives';
  ok $m1.get(0,1) == 12e0, 'element set using postcircumfix lives';
  dies-ok { $m1[3;3] = 12e0 }, 'set element (4,4) using postcircumfix of a (3,3) matrix dies';
  lives-ok { $m1.setall(1) }, 'set all elements';
  is-deeply ($m1.get-row($_) for ^3).flat, (1e0 xx 9), 'all elements set';
  lives-ok { $m1.zero }, 'set all elements to zero';
  is-deeply ($m1.get-row($_) for ^3).flat, (0e0 xx 9), 'all elements zeroed';
  lives-ok { $m1.identity }, 'identity';
  is-deeply ($m1.get-row($_) for ^3).flat, (1e0, 0e0, 0e0, 0e0, 1e0, 0e0, 0e0, 0e0, 1e0), 'identity matrix';
};

subtest 'matrix IO' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  $m1.setall(1);
  my $filename = 't/matrix.dat';
  LEAVE { with $filename.IO { .unlink if .e } }
  lives-ok { $m1.write($filename) }, 'can write a matrix to file';
  lives-ok { $m1.read($filename) }, 'can read a matrix from file';
  lives-ok { $m1.printf($filename, "%3.1lf") }, 'can print a matrix to file';
  is-deeply $filename.IO.lines, '1.0' xx 9, 'matrix file';
};

subtest 'views' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(4));
  $m1.setall(1);
  my Math::Libgsl::Matrix::View $mv1 .= new;
  my Math::Libgsl::Matrix $m2 = $mv1.submatrix($m1, 1, 1, 2, 2);
  $m2.setall(12);
  is-deeply ($m1.get-row($_) for ^3).flat, (1e0, 1e0, 1e0, 1e0, (1e0, 12e0, 12e0, 1e0) xx 2).flat, 'modify matrix using submatrix';
  dies-ok  { $mv1.submatrix($m1, 3, 2, 1, 1) }, 'submatrix 4x4 of a 3x4 matrix dies';
  dies-ok  { $mv1.submatrix($m1, 2, 4, 1, 1) }, 'submatrix 3x5 of a 3x4 matrix dies';
  dies-ok  { $mv1.submatrix($m1, 2, 3, 2, 1) }, 'submatrix 4x4 of a 3x4 matrix dies';
  my Math::Libgsl::Vector $v1 .= new(:size(9));
  $v1.set($_, $_.Num) for ^9;
  $m2 = $mv1.vector($v1, 3, 3);
  is-deeply ($m2.get-row($_) for ^3).flat, (0e0..8e0).flat, 'matrix view of a vector';
  $m2 = $mv1.vector-tda($v1, 2, 2, 3);
  is-deeply ($m2.get-row($_) for ^2).flat, (0e0, 1e0, 3e0, 4e0), 'matrix view with tda from a vector';
  dies-ok { $mv1.vector-tda($v1, 3, 4, 3) }, "can create a matrix view with tda if tda<n2";
  $m1.set(0, $_, $_) for ^4;
  $m1.set(1, 0, 1);
  $m1.set(2, 0, 2);
  my Math::Libgsl::Vector $v2;
  my Math::Libgsl::Vector::View $vv .= new;
  $v2 = $m1.row-view($vv, 0);
  is-deeply $v2[^4], (^4)».Num.flat, 'row view';
  $v2 = $m1.col-view($vv, 0);
  is-deeply $v2[^3], (^3)».Num.flat, 'col view';
  $v2 = $m1.subrow-view($vv, 0, 1, 2);
  is-deeply $v2[^2], (1e0, 2e0), 'subrow view';
  $v2 = $m1.subcol-view($vv, 0, 1, 2);
  is-deeply $v2[^2], (1e0, 2e0), 'subcol view';
  $m1.set(0, 0, 1);
  $v2 = $m1.diagonal-view($vv);
  is-deeply $v2[^3], (1e0, 12e0, 12e0), 'diagonal view';
  $v2 = $m1.subdiagonal-view($vv, 1);
  is-deeply $v2[^2], (1e0, 12e0), 'subdiagonal view';
  dies-ok { $m1.subdiagonal-view($vv, 6) }, "can't create a subdiagonal view with k>min(size1,size2)";
  $v2 = $m1.superdiagonal-view($vv, 1);
  is-deeply $v2[^2], (1e0, 12e0), 'superdiagonal view';
  my @a = 1, 2, 3, 4, 5, 6;
  my $array = prepmat @a;
  $m2 = $mv1.array($array, 2, 3);
  is-deeply ($m2.get-row($_) for ^2).flat, (1e0..6e0).flat, 'matrix view of an array';
  $m2 = $mv1.array-tda($array, 2, 3, 4);
  is-deeply ($m2.get-row($_) for ^2).flat, (1e0, 2e0, 3e0, 5e0, 6e0, 0e0).flat, 'matrix view of an array with tda';
  dies-ok { $mv1.array-tda($array, 2, 3, 2) }, "can't create a view with tda if tda<size2";
  $array = num64-prepmat @a;
  $m2 = $mv1.array($array, 2, 3);
  is-deeply ($m2.get-row($_) for ^2).flat, (1e0..6e0).flat, 'num64-prepmat works';
  my $list := (1, 2, 3, 4, 5, 6);
  my Math::Libgsl::Matrix $matrix .= new: 2, 3;
  array-mat(-> $mat { $matrix.copy($mat); $matrix[0;0] = 0 }, 2, 3, $list);
  is-deeply ($matrix.get-row($_) for ^2).flat, (0e0, 2e0..6e0).flat, 'matrix view of a list as a closure';
  my @array = 1, 2, 3, 4, 5, 6;
  array-mat(sub ($mat) { $matrix.copy($mat); $matrix[0;0] = 0 }, 2, 3, @array);
  is-deeply ($matrix.get-row($_) for ^2).flat, (0e0, 2e0..6e0).flat, 'matrix view of an array as an anonymous sub';
  array-tda-mat(-> $mat { $matrix.copy($mat); $matrix[0;0] = 0 }, 2, 3, 4, @array);
  is-deeply ($matrix.get-row($_) for ^2).flat, (0e0, 2e0, 3e0, 5e0, 6e0, 0e0).flat, 'matrix view of an array with tda';
  array-tda-mat(-> $mat { $matrix.copy($mat); $matrix[0;0] = 0 }, 2, 3, 4, $list);
  is-deeply ($matrix.get-row($_) for ^2).flat, (0e0, 2e0, 3e0, 5e0, 6e0, 0e0).flat, 'matrix view of a list with tda';
};

subtest 'copy & swap' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  $m1.setall(1);
  my Math::Libgsl::Matrix $m2 .= new(:size1(3), :size2(3));
  lives-ok { $m2.copy($m1) }, 'can copy a matrix';
  is-deeply ($m2.get-row($_) for ^3).flat, (1e0 xx 9), 'matrix copy';
  $m1.setall(2);
  lives-ok { $m2.swap($m1) }, 'can swap two matrices';
  is-deeply ($m2.get-row($_) for ^3).flat, (2e0 xx 9), 'matrix swap';
  my Math::Libgsl::Matrix $m3 .= new(:size1(4), :size2(4));
  $m3.setall(1);
  my Math::Libgsl::Matrix $m4 .= new(:size1(4), :size2(4));
  lives-ok { $m4.tricpy($m3, CblasUpper, CblasUnit) }, 'can tricpy';
  is-deeply ($m4.get-row($_) for ^4).flat, (0e0, 1e0, 1e0, 1e0, 0e0, 0e0, 1e0, 1e0, 0e0, 0e0, 0e0, 1e0, 0e0, 0e0, 0e0, 0e0), 'tricpy matrix';
};

subtest 'get & set rows & columns' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(4));
  my Math::Libgsl::Vector $v1 .= new(:size(4));
  $v1.set($_, $_.Num) for ^4;
  lives-ok { $m1.set-row(0, $v1) }, 'can set a row from a vector';
  is-deeply $m1.get-row(0), (^4)».Num.flat.Array, 'set row';
  my @row = ^4;
  lives-ok { $m1.set-row(0, @row) }, 'can set a row from an array';
  is-deeply $m1.get-row(0), (^4)».Num.flat.Array, 'set row';
  my Math::Libgsl::Vector $v2 .= new(:size(3));
  $v2.set($_, $_.Num) for ^3;
  lives-ok { $m1.set-col(0, $v2) }, 'can set a row from a vector';
  is-deeply $m1.get-col(0), (^3)».Num.flat.Array, 'set col';
  my @col = ^3;
  lives-ok { $m1.set-col(0, @col) }, 'can set a col from an array';
  is-deeply $m1.get-col(0), (^3)».Num.flat.Array, 'set col';
};

subtest 'exchange rows & columns' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(4));
  $m1.set-row(0, ^4);
  $m1.set-row(1, 4..7);
  lives-ok { $m1.swap-rows(0, 1) }, 'can swap rows';
  is-deeply $m1.get-row(1), (^4)».Num.flat.Array, 'swap rows';
  lives-ok { $m1.swap-cols(0, 1) }, 'can swap cols';
  is-deeply $m1.get-col(0), [5e0, 1e0, 0e0], 'swap cols';
  dies-ok { $m1.swap-rowcol(0, 1) }, "can't swap row & col when not square";
  my Math::Libgsl::Matrix $m2 .= new(:size1(3), :size2(3));
  $m2.set-row(0, ^3);
  $m2.set-row(1, 3..5);
  $m2.set-row(2, 6..8);
  lives-ok { $m2.swap-rowcol(0, 0) }, 'can swap row & col';
  is-deeply $m2.get-row(0), [0e0, 3e0, 6e0], 'swap row & col';
  my Math::Libgsl::Matrix $m3 .= new(:size1(4), :size2(3));
  lives-ok { $m3.copy-transpose($m1) }, 'can copy and transpose';
  is-deeply ($m3.get-row($_) for ^4).flat, (5e0, 1e0, 0e0, 4e0, 0e0, 0e0, 6e0, 2e0, 0e0, 7e0, 3e0, 0e0), 'copy and transpose';
  lives-ok { $m2.transpose() }, 'can transpose';
  is-deeply ($m2.get-row($_) for ^3).flat, (0e0, 1e0, 2e0, 3e0, 4e0, 7e0, 6e0, 5e0, 8e0), 'transpose';
  my Math::Libgsl::Matrix $m4 .= new(:size1(4), :size2(4));
  $m4.setall(1);
  my Math::Libgsl::Matrix $m5 .= new(:size1(4), :size2(4));
  lives-ok { $m5.transpose-tricpy($m4, CblasUpper, CblasUnit) }, 'can tricpy and transpose';
  is-deeply ($m5.get-row($_) for ^4).flat, (0e0, 0e0, 0e0, 0e0, 1e0, 0e0, 0e0, 0e0, 1e0, 1e0, 0e0, 0e0, 1e0, 1e0, 1e0, 0e0), 'tricpy and transpose matrix';
};

subtest 'matrix operations' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  my Math::Libgsl::Matrix $m2 .= new(:size1(3), :size2(3));
  $m1.set-row($_, ^3) for ^3;
  $m2.set-row($_, ^3) for ^3;
  $m1.add($m2);
  is-deeply ($m1.get-row($_) for ^3).flat, ((0e0, 2e0, 4e0) xx 3).flat, 'can add';
  $m1.sub($m2);
  is-deeply ($m1.get-row($_) for ^3).flat, ((0e0, 1e0, 2e0) xx 3).flat, 'can sub';
  $m1.mul($m2);
  is-deeply ($m1.get-row($_) for ^3).flat, ((0e0, 1e0, 4e0) xx 3).flat, 'can mul';
  $m1.set($_, 0, 1) for ^3;
  $m2.set($_, 0, 1) for ^3;
  $m1.div($m2);
  is-deeply ($m1.get-row($_) for ^3).flat, ((1e0, 1e0, 2e0) xx 3).flat, 'can div';
  $m1.scale(2);
  is-deeply ($m1.get-row($_) for ^3).flat, ((2e0, 2e0, 4e0) xx 3).flat, 'can scale';
  $m1.add-constant(1);
  is-deeply ($m1.get-row($_) for ^3).flat, ((3e0, 3e0, 5e0) xx 3).flat, 'can add a constant';
  if $gsl-version > v2.6 {
    my Math::Libgsl::Vector $v1 .= new(:size(3));
    $v1.set($_, $_.Num) for ^3;
    $m1.scale-rows($v1);
    is-deeply ($m1.get-row($_) for ^3).flat, (0e0, 0e0, 0e0, 3e0, 3e0, 5e0, 6e0, 6e0, 10e0).flat, 'can scale-rows';
    $m1.set-row($_, ^3) for ^3;
    $m1.scale-columns($v1);
    is-deeply ($m1.get-row($_) for ^3).flat, (0e0, 1e0, 4e0, 0e0, 1e0, 4e0, 0e0, 1e0, 4e0).flat, 'can scale-columns';
  }
};

subtest 'matrix min & max' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  $m1.set-row(0, ^3);
  $m1.set-row(1, 3..5);
  $m1.set-row(2, 6..8);
  ok $m1.max == 8e0, 'max element';
  ok $m1.min == 0e0, 'min element';
  is-deeply $m1.minmax, (0e0, 8e0), 'min & max elements';
  is-deeply $m1.max-index, (2, 2), 'index of max element';
  is-deeply $m1.min-index, (0, 0), 'index of min element';
  is-deeply $m1.minmax-index, (0, 0, 2, 2), 'index of min & max elements';
};

subtest 'matrix properties' => {
  my Math::Libgsl::Matrix $m1 .= new(:size1(3), :size2(3));
  ok $m1.is-null, 'all elements are 0';
  $m1.setall(2e0);
  ok $m1.is-pos, 'all elements are >0';
  $m1.set(0, 0, 0e0);
  nok $m1.is-pos, 'fails if not all elements are >0';
  $m1.setall(-2e0);
  ok $m1.is-neg, 'all elements are <0';
  my Math::Libgsl::Matrix $m2 .= new(:size1(3), :size2(3));
  $m2.setall(-2e0);
  ok $m1.is-equal($m2), 'two matrices are equal';
  if $gsl-version > v2.6 {
    $m1.setall(2);
    ok $m1.norm1 == 6, 'matrix 1-norm';
  }
};

subtest 'matrix to list' => {
  my Math::Libgsl::Matrix $m .= new(:size1(3), :size2(3)).identity;
  is-deeply $m.list».Int, ((1, 0, 0), (0, 1, 0), (0, 0, 1)), 'list method coerces to a List of Lists';
  my @l = $m.list;
  is-deeply @l».Int, [(1, 0, 0), (0, 1, 0), (0, 0, 1)], 'assign a matrix to a Raku array';
};

subtest 'gist & put' => {
  my Math::Libgsl::Matrix $m .= new(:size1(200), :size2(200)).identity;

  for ^200 X ^200 -> ($i, $j) {
    $m[$i;$j] = $j + $i * 200;
  }

  ok $m.gist eq '(' ~ do for ^10 -> $i { '(' ~ ((^10) »+» ($i * 200)).join(' ') ~ " ...)\n" } ~ "...)",
    'gist returns the first 10 elements of the first 10 rows';

  sub capture_stdout(Callable $code) {
    my $result;
    my $*OUT = class {
      method print(*@args) {
        $result ~= @args.join;
      }
      method flush {}
    }
    $code.();
    $result;
  }

  my $out = capture_stdout({ $m.put });
  ok $out eq (^40000).join(' ') ~ "\n", 'method put';
  $out = capture_stdout({ put $m });
  ok $out eq (^40000).join(' ') ~ "\n", 'sub put';
};

done-testing;


