### ----------------------------------------------------
### -- Math::Matrix
### -- Licenses: Artistic-2.0
### -- Authors: Pierre Vigier, Herbert Breunung
### -- File: lib/Math/Matrix.rakumod
### ----------------------------------------------------
### -- Chunk 2 of 2
i][$_] * @L[$k][$_] for 0 ..^ $k ;
            @L[$i][$k] /= @L[$k][$k];
        }
    }
    for ^$!row-count X ^$!column-count -> ($r, $c) { @L[$r][$c] = 0 if $r < $c }
    return Math::Matrix.new-lower-triangular( @L ) unless $diagonal;
    my @D;
    for ^$!row-count -> $r {
        @D[$r] = @L[$r][$r] ** 2 ;
        @L[$r][$r] = 1;
    }
    Math::Matrix.new-lower-triangular( @L ), Math::Matrix.new-diagonal( @D );
}

################################################################################
# end of decompositions - start matrix math operations
################################################################################


method equal(Math::Matrix:D: Math::Matrix $b --> Bool)           { @!rows ~~ $b!rows }
multi method ACCEPTS(Math::Matrix:D: Math::Matrix:D $b --> Bool) { self.equal( $b )  }

# add matrix
multi method add(Math::Matrix:D: Str $b --> Math::Matrix:D ) { self.add( Math::Matrix.new( $b ) ) }
multi method add(Math::Matrix:D:     @b --> Math::Matrix:D ) { self.add( Math::Matrix.new( @b ) ) }
multi method add(Math::Matrix:D: Math::Matrix $b where {self.size eqv $b.size} --> Math::Matrix:D ) {
    my @sum;
    for ^$!row-count X ^$!column-count -> ($r, $c) {
        @sum[$r][$c] = @!rows[$r][$c] + $b!rows[$r][$c];
    }
    Math::Matrix.new( @sum );
}
# add vector
multi method add(Math::Matrix:D: @v where {@v.all ~~ Numeric}, Int :$row! --> Math::Matrix:D) {
    fail "Matrix has $!column-count columns, but got "~ +@v ~ "element row." unless $!column-count == +@v;
    self!check-row-index($row);
    my @m = self!clone-cells;
    @m[$row] = @m[$row] <<+>> @v;
    Math::Matrix.new( @m );
}
multi method add(Math::Matrix:D: @v where {@v.all ~~ Numeric}, Int :$column! --> Math::Matrix:D ) {
    self!check-column-index($column);
    fail "Matrix has $!row-count rows, but got "~ +@v ~ "element column." unless $!row-count == +@v;
    my @m = self!clone-cells;
    @v.keys.map:{ @m[$_][$column] += @v[$_] };
    Math::Matrix.new( @m );
}
# add scalar
multi method add(Math::Matrix:D: Numeric $s, --> Math::Matrix:D )                         { self.map( *  + $s ) }
multi method add(Math::Matrix:D: Numeric $s, Int :$row, Int :$column --> Math::Matrix:D ) {
    self!check-row-index($row)       if $row.defined;
    self!check-column-index($column) if $column.defined;
    my @sum = self!clone-cells;
    for   ($row.defined ?? ($row..$row) !! ^$!row-count)
        X ($column.defined ?? ($column..$column) !! ^$!column-count)
       -> ($r, $c) { @sum[$r][$c] += $s }
    Math::Matrix.new( @sum );
}

# multiply matrix
multi method multiply(Math::Matrix:D: Str $b --> Math::Matrix:D ) { self.multiply( Math::Matrix.new( $b ) ) }
multi method multiply(Math::Matrix:D:     @b --> Math::Matrix:D ) { self.multiply( Math::Matrix.new( @b ) ) }
multi method multiply(Math::Matrix:D: Math::Matrix $b where {self.size eqv $b.size} --> Math::Matrix:D ) {
    my @product;
    for ^$!row-count X ^$!column-count -> ($r, $c) { @product[$r][$c] = @!rows[$r][$c] * $b!rows[$r][$c] }
    Math::Matrix.new( @product );
}
# multiply vector
multi method multiply(Math::Matrix:D: @v where {@v.all ~~ Numeric}, Int :$row! --> Math::Matrix:D) {
    fail "Matrix has $!column-count columns, but got "~ +@v ~ "element row." unless $!column-count == +@v;
    self!check-row-index($row);
    my @m = self!clone-cells;
    @m[$row] = @m[$row] <<*>> @v;
    Math::Matrix.new( @m );
}
multi method multiply(Math::Matrix:D: @v where {@v.all ~~ Numeric}, Int :$column! --> Math::Matrix:D ) {
    self!check-column-index($column);
    fail "Matrix has $!row-count rows, but got "~ +@v ~ "element column." unless $!row-count == +@v;
    my @m = self!clone-cells;
    @v.keys.map:{ @m[$_][$column] *= @v[$_] };
    Math::Matrix.new( @m );
}
# multiply scalar
multi method multiply(Math::Matrix:D: Numeric $f --> Math::Matrix:D )                          { self.map( *  * $f  ) }
multi method multiply(Math::Matrix:D: Numeric $f, Int :$row, Int :$column --> Math::Matrix:D ) {
    self!check-row-index($row)       if $row.defined;
    self!check-column-index($column) if $column.defined;
    my @product = self!clone-cells;
    for   ($row.defined ?? ($row..$row) !! ^$!row-count)
        X ($column.defined ?? ($column..$column) !! ^$!column-count)
       -> ($r, $c) { @product[$r][$c] *= $f }
    Math::Matrix.new( @product );
}

method dot-product(Math::Matrix:D: Math::Matrix $b --> Math::Matrix:D ) {
    fail "Number of columns of the second matrix is different from number of rows of the first operand"
        unless $!column-count == $b!row-count;
    my @product;
    for ^$!row-count X ^$b!column-count -> ($r, $c) {
        @product[$r][$c] += @!rows[$r][$_] * $b!rows[$_][$c] for ^$b!row-count;
    }
    Math::Matrix.new( @product );
}

method tensor-product(Math::Matrix:D: Math::Matrix $b  --> Math::Matrix:D) {
    my @product;
    for @!rows -> $arow {
        for $b!rows -> $brow {
            @product.push([ ($arow.list.map: { $brow.flat >>*>> $_ }).flat ]);
        }
    }
    Math::Matrix.new( @product );
}

################################################################################
# end of matrix math operations - start list like operations
################################################################################

method elems (Math::Matrix:D: --> Int)           {  $!row-count * $!column-count }

method elem (Math::Matrix:D: Range $r --> Bool)  {         # is every element value element in the set/range
    self.list.map: {return False unless $_ ~~ $r};
    True;
}

multi method cont (Math::Matrix:D: Numeric $e  --> Bool) { # matrix contains element ?
    self.list.map: {return True if $_ == $e};
    False;
}
multi method cont (Math::Matrix:D: Range $r  --> Bool) {   # is any element value in this set/range
    self.list.map: {return True if $_ ~~ $r};
    False;
}

method map(Math::Matrix:D: &coderef, Range :$rows = ^$!row-count,
                                     Range :$columns = ^$!column-count --> Math::Matrix:D) {
    self!check-index($rows.min, $columns.min);
    self!check-row-index($rows.minmax[1]) unless $rows.max          == Inf;
    self!check-column-index($columns.minmax[1]) unless $columns.max == Inf;
    my @r = $rows.max    == Inf ?? ($rows.min    .. $!row-count-1).list    !! $rows.list;
    my @c = $columns.max == Inf ?? ($columns.min .. $!column-count-1).list !! $columns.list;
    my @m;
    for @r X @c                        -> ($r, $c) { @m[$r][$c] = &coderef(@!rows[$r][$c])}
    for ^$!row-count X ^$!column-count -> ($r, $c) { @m[$r][$c] //= @!rows[$r][$c] }
    Math::Matrix.new( @m );
}

method map-with-index(Math::Matrix:D: &coderef, Range :$rows = ^$!row-count,
                                                Range :$columns = ^$!column-count --> Math::Matrix:D) {
    fail "block has to receive between one and three arguments" unless &coderef.arity ~~ 1..3;
    self!check-index($rows.min, $columns.min);
    self!check-row-index($rows.minmax[1]) unless $rows.max          == Inf;
    self!check-column-index($columns.minmax[1]) unless $columns.max == Inf;
    my @r = $rows.max    == Inf ?? ($rows.min    .. $!row-count-1).list    !! $rows.list;
    my @c = $columns.max == Inf ?? ($columns.min .. $!column-count-1).list !! $columns.list;
    my @m;
    if    &coderef.arity == 1 {for @r X @c -> ($r, $c) { @m[$r][$c] = &coderef($r) }}
    elsif &coderef.arity == 2 {for @r X @c -> ($r, $c) { @m[$r][$c] = &coderef($r, $c) }}
    elsif &coderef.arity == 3 {for @r X @c -> ($r, $c) { @m[$r][$c] = &coderef($r, $c, @!rows[$r][$c]) }}
    for ^$!row-count X ^$!column-count     -> ($r, $c) { @m[$r][$c] //= @!rows[$r][$c] }
    Math::Matrix.new( @m );
}


method reduce(        Math::Matrix:D: &coderef) {(@!rows.map: {$_.flat}).flat.reduce( &coderef )}
method reduce-rows   (Math::Matrix:D: &coderef) { @!rows.map: { $_.flat.reduce( &coderef) }}
method reduce-columns(Math::Matrix:D: &coderef) {(^$!column-count).map: { self.column($_).reduce( &coderef )}}

################################################################################
# end of list like operations - start structural matrix operations
################################################################################

multi method move-row (Math::Matrix:D: Pair $p --> Math::Matrix:D) {
    self.move-row($p.key, $p.value)
}
multi method move-row (Math::Matrix:D: Int $from, Int $to --> Math::Matrix:D) {
    self!check-row-indices([$from, $to]);
    return self if $from == $to;
    my @rows = (^$!row-count).list;
    @rows.splice($to,0,@rows.splice($from,1));
    self.submatrix( rows => @rows, columns => (^$!column-count).list);
}

multi method move-column (Math::Matrix:D: Pair $p --> Math::Matrix:D) {
    self.move-column($p.key, $p.value)
}
multi method move-column (Math::Matrix:D: Int $from, Int $to --> Math::Matrix:D) {
    self!check-column-indices([$from, $to]);
    return self if $from == $to;
    my @cols = (^$!column-count).list;
    @cols.splice($to,0,@cols.splice($from,1));
    self.submatrix( rows => (^$!row-count).list, columns => @cols);
}

method swap-rows (Math::Matrix:D: Int $rowa, Int $rowb --> Math::Matrix:D) {
    self!check-row-indices([$rowa, $rowb]);
    return self if $rowa == $rowb;
    my @rows = (^$!row-count).list;
    (@rows.[$rowa], @rows.[$rowb]) = (@rows.[$rowb], @rows.[$rowa]);
    self.submatrix( rows => @rows, columns => (^$!column-count).list);
}

method swap-columns (Math::Matrix:D: Int $cola, Int $colb --> Math::Matrix:D) {
    self!check-column-indices([$cola, $colb]);
    return self if $cola == $colb;
    my @cols = (^$!column-count).list;
    (@cols.[$cola], @cols.[$colb]) = (@cols.[$colb], @cols.[$cola]);
    self.submatrix( rows => (^$!row-count).list, columns => @cols);
}

multi method splice-rows(Math::Matrix:D: Int $row, Int $elems, Math::Matrix $replacement --> Math::Matrix:D){
    self.splice-rows($row, $elems, $replacement.Array );
}
multi method splice-rows(Math::Matrix:D: Int $row, Int $elems = ($!row-count - $row), Array $replacement = [] --> Math::Matrix:D){
    my $pos = $row >= 0 ?? $row !! $!row-count + $row + 1;
    fail "Row index (first parameter) is outside of matrix size!" unless 0 <= $pos <= $!row-count;
    fail "Number of elements to delete (second parameter) has to be zero or more!)" if $elems < 0;
    if $replacement.elems > 0 {
        fail "Number of columns in and original matrix and replacement has to be same" unless $replacement[0].elems == $!column-count;
        Math::Matrix::ArrayOfArray::check-data( @$replacement );
    }
    my @m = self!clone-cells;
    @m.splice($pos, $elems, $replacement.list);
    Math::Matrix.new(@m);
}


multi method splice-columns(Math::Matrix:D: Int $col, Int $elems, Math::Matrix $replacement --> Math::Matrix:D){
    self.splice-columns($col, $elems, $replacement.Array );
}
multi method splice-columns(Math::Matrix:D: Int $col, Int $elems = ($!column-count - $col), Array $replacement = ([[] xx $!row-count]) --> Math::Matrix:D){
    my $pos = $col >= 0 ?? $col !! $!column-count + $col + 1;
    fail "Column index (first parameter) is outside of matrix size!" unless 0 <= $pos <= $!column-count;
    fail "Number of elements to delete (second parameter) has to be zero or more!)" if $elems < 0;
    fail "Number of rows in original matrix and replacement has to be same" unless $replacement.elems == $!row-count;
    Math::Matrix::ArrayOfArray::check-data( @$replacement );
    my @m = self!clone-cells;
    @m.keys.map:{ @m[$_].splice($pos, $elems, $replacement[$_]) };
    Math::Matrix.new(@m);
}

################################################################################
# end of structural matrix operations - start operators
################################################################################

multi sub prefix:<@>( Math::Matrix:D $m --> Array)          is export { $m.Array }
multi sub prefix:<%>( Math::Matrix:D $m --> Hash)           is export { $m.Hash }
multi sub prefix:<->( Math::Matrix:D $m --> Math::Matrix:D) is export { $m.negated }

multi sub circumfix:<｜ ｜>( Math::Matrix:D $m --> Numeric) is equiv(&prefix:<!>) is export { $m.determinant }
multi sub circumfix:<‖ ‖>( Math::Matrix:D $m --> Numeric)  is equiv(&prefix:<!>) is export { $m.norm }


multi sub infix:<+>( Math::Matrix:D $a, Numeric $n        --> Math::Matrix:D) is export { $a.add($n) }
multi sub infix:<+>( Numeric $n,        Math::Matrix:D $a --> Math::Matrix:D) is export { $a.add($n) }
multi sub infix:<+>( Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is export { $a.add($b) }

multi sub infix:<->( Numeric $n,        Math::Matrix:D $a --> Math::Matrix:D) is export { $a.negated.add($n) }
multi sub infix:<->( Math::Matrix:D $a, Numeric $n        --> Math::Matrix:D) is export { $a.add(-$n)  }
multi sub infix:<->( Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is export { $a.subtract($b) }

multi sub infix:<*>( Math::Matrix:D $a, Numeric $n        --> Math::Matrix:D) is export { $a.multiply($n) }
multi sub infix:<*>( Numeric $n,        Math::Matrix:D $a --> Math::Matrix:D) is export { $a.multiply($n) }
multi sub infix:<*>( Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is export { $a.multiply($b) }
multi sub infix:<**>( Math::Matrix:D $a where { $a.is-square }, Int $e --> Math::Matrix:D) is export {
    return Math::Matrix.new-identity( $a!row-count ) if $e == 0;
    my $p = $a.clone;
    $p = $p.dot-product( $a ) for 2 .. abs $e;
    $p = $p.inverted          if  $e < 0;
    $p;
}

multi sub infix:<⋅>(  Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is tighter(&infix:<*>) is export { $a.dot-product( $b ) }
multi sub infix:<dot>(Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is equiv(&infix:<⋅>)   is export { $a.dot-product( $b ) }
multi sub infix:<÷>(  Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is equiv(&infix:<⋅>)   is export { $a.dot-product( $b.inverted ) }

multi sub infix:<⊗>( Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is equiv(&infix:<x>) is export { $a.tensor-product( $b ) }
multi sub infix:<X*>( Math::Matrix:D $a, Math::Matrix:D $b --> Math::Matrix:D) is equiv(&infix:<x>) is export { $a.tensor-product( $b ) }

multi sub prefix:<MM>(Str   $m --> Math::Matrix:D) is tighter(&postcircumfix:<[ ]>) is export(:MM) { Math::Matrix.new($m) }
multi sub prefix:<MM>(List  $m --> Math::Matrix:D) is tighter(&postcircumfix:<[ ]>) is export(:MM) { Math::Matrix.new(@$m) }

# vim: expandtab shiftwidth=4


