# -*- tcl -*-
# linalg.test --
#    Tests for the linear algebra package
#
# NOTE:
#    Comparison by numbers, not strings, needed!
#
# TODO:
#    Tests for:
#    - show, angle
#    - solveGaussBand, solveTriangularBand
#    - mkHilbert and so on
#    - matmul

# -------------------------------------------------------------------------

source [file join \
	[file dirname [file dirname [file join [pwd] [info script]]]] \
	devtools testutilities.tcl]

testsNeedTcl     8.4
testsNeedTcltest 2.1

support {
    useLocal math.tcl math
}
testing {
    useLocal linalg.tcl math::linearalgebra
}

# -------------------------------------------------------------------------

namespace import -force ::math::linearalgebra::*

set old_precision $::tcl_precision
set ::tcl_precision 17

#
# Matching procedure - flatten the lists
#
proc matchNumbers {expected actual} {
    set match 1
    if { [llength [lindex $expected 0]] > 1 } {
        foreach a $actual e $expected {
            set match [matchNumbers $e $a]
            if { $match == 0 } {
                break
            }
        }
    } else {
        foreach a $actual e $expected {
            if {abs($a-$e) > 0.1e-6} {
                set match 0
                break
            }
        }
    }
    return $match
}

customMatch numbers matchNumbers

test dimshape-1.0 "dimension of scalar" -body {
    dim 1
} -result 0

test dimshape-1.1 "dimension of vector" -body {
    dim {1 2 3}
} -result 1

test dimshape-1.2 "dimension of matrix" -body {
    dim { {1 2 3} {4 5 6} }
} -result 2

test dimshape-2.0 "shape of scalar" -body {
    shape 1
} -result {}

test dimshape-2.1 "shape of vector" -body {
    shape {1 2 3}
} -result 3

test dimshape-2.2 "shape of matrix" -body {
    shape { {1 2 3} {4 5 6} }
} -result {2 3}

test symmetric-1.0 "non-symmetric matrix" -body {
    symmetric { {1 2 3} {4 5 6} {7 8 9}}
} -result 0

test symmetric-1.1 "symmetric matrix" -body {
    symmetric { {1 2 3} {2 1 4} {3 4 1}}
} -result 1

test symmetric-1.2 "non-square matrix" -body {
    symmetric { {1 2 3} {2 1 4}}
} -result 0

test norm-1.0 "one-norm - 5 components" -match numbers -body {
    norm {1 2 3 0 -1} 1
} -result 7.0

test norm-1.1 "one-norm - 2 components" -match numbers -body {
    norm {1 -1} 1
} -result 2.0

test norm-1.2 "two-norm - 5 components" -match numbers -body {
    norm {1 2 3 0 -1} 2
} -result [expr {sqrt(15)}]

test norm-1.3 "two-norm - 2 components" -match numbers -body {
    norm {1 -1} 2
} -result [expr {sqrt(2)}]

test norm-1.4 "two-norm - no underflow" -match numbers -body {
    norm {3.0e-140 -4.0e-140} 2
} -result 5.0e-140

test norm-1.5 "two-norm - no overflow" -match numbers -body {
    norm {3.0e140 -4.0e140} 2
} -result 5.0e140

test norm-1.6 "max-norm - 5 components" -match numbers -body {
    norm {1 2 3 0 -4} max
} -result 4

test norm-1.7 "max-norm - 2 components" -match numbers -body {
    norm {1 -1} max
} -result 1

test norm-2.0 "matrix-norm - 2x2 - max" -match numbers -body {
    normMatrix {{1 -1} {1 1}} max
} -result 1

test norm-2.1 "matrix-norm - 2x2 - 1" -match numbers -body {
    normMatrix {{1 -1} {1 1}} 1
} -result 4

test norm-2.2 "matrix-norm - 2x2 - 2" -match numbers -body {
    normMatrix {{1 -1} {1 1}} 2
} -result 2

test norm-3.0 "statistical normalisation - vector" -match numbers -body {
    normalizeStat {1 0 0 0}
} -result {1.5 -0.5 -0.5 -0.5}

test norm-3.1 "statistical normalisation - matrix" -match numbers -body {
    normalizeStat {{1 0 0 0} {0 0 0 1} {0 1 1 0} {0 0 0 0}}
} -result {{ 1.5 -0.5 -0.5 -0.5}
           {-0.5 -0.5 -0.5  1.5}
           {-0.5  1.5  1.5 -0.5}
           {-0.5 -0.5 -0.5 -0.5}}

test dotproduct-1.0" "dot-product - 2 components" -match numbers -body {
    dotproduct {1 -1} {1 -1}
} -result 2.0

test dotproduct-1.1" "dot-product - 5 components" -match numbers -body {
    dotproduct {1 2 3 4 5} {5 4 3 2 1}
} -result [expr {5.0+8+9+8+5}]

test unitlength-1.0" "unitlength - 2 components" -match numbers -body {
    unitLengthVector {3 4}
} -result {0.6 0.8}

test unitlength-1.1" "unitlength - 4 components" -match numbers -body {
    unitLengthVector {1 1 1 1}
} -result {0.5 0.5 0.5 0.5}

test axpy-1.0 "axpy - vectors" -body {
    axpy 2 {1 -1} {2 -2}
} -result {4 -4}

test axpy-1.1 "axpy - matrices" -body {
    axpy 2 { {1 -1} {2 -2} {3 4} {-3 4} } \
           { {5 -5} {5 -5} {6 6} {-6 6} }
} -result {{7 -7} {9 -9} {12 14} {-12 14}}

test add-1.0 "add - vectors" -body {
    add {1 -1} {2 -2}
} -result {3 -3}

test add-1.1 "add - matrices" -body {
    add { {1 -1} {2 -2} {3 4} {-3 4} } \
        { {5 -5} {5 -5} {6 6} {-6 6} }
} -result {{6 -6} {7 -7} {9 10} {-9 10}}

test sub-1.0 "sub - vectors" -body {
    sub {1 -1} {2 -2}
} -result {-1 1}

test sub-1.1 "sub - matrices" -body {
    sub { {1 -1} {2 -2} {3 4} {-3 4} } \
        { {5 -5} {5 -5} {6 6} {-6 6} }
} -result {{-4 4} {-3 3} {-3 -2} {3 -2}}

test scale-1.0 "scale - vectors" -body {
    scale 3 {2 -2}
} -result {6 -6}

test scale-1.1 "scale - matrices" -body {
    scale 3 { {5 -5} {5 -5} {6 6} {-6 6} }
} -result {{15 -15} {15 -15} {18 18} {-18 18}}

test make-1.0 "mkVector - create a null vector" -body {
    mkVector 3
} -result {0.0 0.0 0.0}

test make-1.1 "mkVector - create a vector with values 1" -body {
    mkVector 3 1.0
} -result {1.0 1.0 1.0}

test make-2.0 "mkMatrix - create a matrix with 3 rows, 2 columns" -body {
    mkMatrix 3 2 2.0
} -result {{2.0 2.0} {2.0 2.0} {2.0 2.0}}

test make-2.1 "mkMatrix - create a matrix with 2 rows, 3 columns" -body {
    mkMatrix 2 3 1.0
} -result {{1.0 1.0 1.0} {1.0 1.0 1.0}}

test make-3.0 "mkIdentity - create an identity matrix 2x2" -body {
    mkIdentity 2
} -result {{1.0 0.0} {0.0 1.0}}

test make-3.1 "mkIdentity - create an identity matrix 3x3" -body {
    mkIdentity 3
} -result {{1.0 0.0 0.0} {0.0 1.0 0.0} {0.0 0.0 1.0}}

test make-4.0 "mkDiagonal - create a diagonal matrix 2x2" -body {
    mkDiagonal {2.0 3.0}
} -result {{2.0 0.0} {0.0 3.0}}

test make-4.1 "mkDiagonal - create a diagonal matrix 3x3" -body {
    mkDiagonal {2.0 3.0 4.0}
} -result {{2.0 0.0 0.0} {0.0 3.0 0.0} {0.0 0.0 4.0}}

test getset-1.0 "getrow - get first row from a matrix" -body {
    getrow {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 0
} -result {1 2 3}

test getset-1.1 "getrow - get last row from a matrix" -body {
    getrow {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 3
} -result {10 11 12}

test getset-1.2 "getcol - get first column from a matrix" -body {
    getcol {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 0
} -result {1 4 7 10}

test getset-1.3 "getcol - get last column from a matrix" -body {
    getcol {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 2
} -result {3 6 9 12}

test getset-2.0 "setrow - set first row in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setrow M 0 {3 2 1}
} -result {{3 2 1} {4 5 6} {7 8 9} {10 11 12}}

test getset-2.1 "setrow - set last row in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setrow M 3 {3 2 1}
} -result {{1 2 3} {4 5 6} {7 8 9} {3 2 1}}

test getset-2.2 "setcol - set first column in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setcol M 0 {3 2 1 0}
} -result {{3 2 3} {2 5 6} {1 8 9} {0 11 12}}

test getset-2.3 "setcol - set last column in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setcol M 2 {3 2 1 0}
} -result {{1 2 3} {4 5 2} {7 8 1} {10 11 0}}

test getset-3.0 "getelem - get element (0,0) in a matrix" -body {
    getelem {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 0 0
} -result 1

test getset-3.1 "getelem - set element (1,2) in a matrix" -body {
    getelem {{1 2 3} {4 5 6} {7 8 9} {10 11 12}} 1 2
} -result 6

test getset-3.2 "setelem - set element (0,0) in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setelem M 0 0 100
} -result {{100 2 3} {4 5 6} {7 8 9} {10 11 12}}

test getset-3.3 "setelem - set element (1,2) in a matrix" -body {
    set M {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    setelem M 1 2 100
} -result {{1 2 3} {4 5 100} {7 8 9} {10 11 12}}

test getset-4.0 "getelem - get element 1 from a vector" -body {
    set V {1 2 3}
    getelem $V 1
} -result 2

test getset-4.1 "setelem - set element 1 in a vector" -body {
    set V {1 2 3}
    setelem V 1 4
} -result {1 4 3}

test rotate-1.0 "rotate - over 90 degrees" -body {
    set v1 {1 2 3}
    set v2 {4 5 6}
    rotate 0 1 $v1 $v2
} -result {{-4 -5 -6} {1 2 3}}

test rotate-1.1 "rotate - over 180 degrees" -body {
    set v1 {1 2 3 4 5 6}
    set v2 {7 8 9 10 11 12}
    rotate -1 0 $v1 $v2
} -result {{-1 -2 -3 -4 -5 -6} {-7 -8 -9 -10 -11 -12}}

test matmul-1.0 "multiply matrix - vector" -match numbers -body {
    set v1 {1 2 3}
    set m  {{0 0 1} {0 5 0} {-1 0 0}}
    matmul $m $v1
} -result {3 10 -1}

test matmul-1.1 "multiply vector - matrix" -match numbers -body {
    set v1 {{1 2 3}} ;# Row vector
    set m  {{0 0 1} {0 5 0} {-1 0 0}}
    matmul $v1 $m
} -result {{-3 10 1}}

test matmul-1.2 "multiply matrix - matrix" -match numbers -body {
    set m1 {{0 0 1} {0 5 0} {-1 0 0}}
    set m2 {{0 0 1} {1 5 1} {-1 0 0}}
    matmul $m1 $m2
} -result {{-1 0 0} {5 25 5} {0 0 -1}}

test matmul-1.3 "multiply vector - vector" -match numbers -body {
    set v1 {1 2 3}
    set v2 {4 5 6}
    matmul $v1 $v2
} -result {{4 5 6} {8 10 12} {12 15 18}}

test matmul-1.4 "multiply row vector - column vector" -match numbers -body {
    set v1 [transpose {1 2 3}]
    set v2 {4 5 6}
    matmul $v1 $v2
} -result 32

test matmul-1.5 "multiply column vector - row vector" -match numbers -body {
    set v1 {1 2 3}
    set v2 [transpose {4 5 6}]
    matmul $v1 $v2
} -result {{4 5 6} {8 10 12} {12 15 18}}

test solve-1.1 "solveGauss - 2x2 matrix" -match numbers -body {
    set b {{2 3} {-2 3}}
    set M {{2 3} {-2 3}}
    solveGauss $M $b
} -result {{1 0} {0 1}}

test solve-1.2 "solveGauss - 3x3 matrix" -match numbers -body {
    set b {{2 3 4} {-2 3 4} {1 1 1}}
    set M {{2 3 4} {-2 3 4} {1 1 1}}
    solveGauss $M $b
} -result {{1 0 0} {0 1 0} {0 0 1}}

test solve-1.3 "solveGauss - 3x3 matrix - less trivial" -match numbers -body {
    set b {{6 -3 6} {2 -3 2} {2 -1 2}}
    set M {{2 3 4} {-2 3 4} {1 1 1}}
    solveGauss $M $b
} -result {{1 0 1} {0 -1 0} {1 0 1}}

test orthon-1.0 "orthonormalize columns - 3x3" -match numbers -body {
    set M {{1 1 1}
           {0 1 1}
           {0 0 1}}
    orthonormalizeColumns $M
} -result {{1 0 0}
           {0 1 0}
           {0 0 1}}

test orthon-1.1 "orthonormalize rows - 3x3" -match numbers -body {
    set M {{1 0 0}
           {1 1 0}
           {1 1 1}}
    orthonormalizeRows $M
} -result {{1 0 0}
           {0 1 0}
           {0 0 1}}

test orthon-1.2 "orthonormalize rows - 3x4" -match numbers -body {
    set M {{1 0 0 0}
           {1 1 0 0}
           {1 1 1 0}}
    orthonormalizeRows $M
} -result {{1 0 0 0}
           {0 1 0 0}
           {0 0 1 0}}

#
# The results from the original LA package have been used
# as a benchmark:
#
#
test svd-1.0 "singular value decomposition - 2x2" -match numbers -body {
    set M {{1.0 2.0} {2.0 1.0}}
    determineSVD $M
} -result {
{{0.70710678118654757  0.70710678118654746}
 {0.70710678118654746 -0.70710678118654757}}
 {3.0 1.0}
{{0.70710678118654757 -0.70710678118654746}
 {0.70710678118654746  0.70710678118654757}}
}

test svd-1.1 "singular value decomposition - 10x10" -match numbers -body {
    set M [mkDingdong 10]
    show [lindex [determineSVD $M] 1] %6.4f
} -result {1.5708 1.5708 1.5708 1.5708 1.5708 1.5707 1.5695 1.5521 1.3935 0.6505}




test LA-1.0 "to_LA - vector" -match numbers -body {
    set vector {1 2 3}
    to_LA $vector
} -result {2 3 0 1 2 3}

test LA-1.1 "to_LA - matrix" -match numbers -body {
    set matrix {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}
    to_LA $matrix
} -result {2 4 3 1 2 3 4 5 6 7 8 9 10 11 12}

test LA-2.0 "from_LA - vector" -match numbers -body {
    set vector {2 3 0 1 2 3}
    from_LA $vector
} -result {1 2 3}

test LA-2.1 "from_LA - matrix" -match numbers -body {
    set matrix {2 4 3 1 2 3 4 5 6 7 8 9 10 11 12}
    from_LA $matrix
} -result {{1 2 3} {4 5 6} {7 8 9} {10 11 12}}

test choleski-1.0 "choleski decomposition of Moler matrix" -match numbers -body {
    set matrix [mkMoler 5]
    choleski $matrix
} -result {{1 0 0 0 0} {-1 1 0 0 0} {-1 -1 1 0 0} {-1 -1 -1 1 0} {-1 -1 -1 -1 1}}

test leastsquares-1.0 "Least-squares solution" -match numbers -body {
    #
    # Known relation: z = 1.0 + x + 0.1*y
    # Model this as: z = z0 + x + 0.1*y
    # (The column of 1s allows us to use a non-zero intercept)
    #
    #          z0   x   y     z
    set Ab { {  1  1.0  1.0} 2.1
             {  1  2.0  1.0} 3.1
             {  1  2.0  2.0} 3.2
             {  1  4.0  2.0} 5.2
             {  1  4.0 22.0} 7.2
             {  1  5.0 -2.0} 5.8 }

    set A {}
    set b {}
    foreach {Ar br} $Ab {
        lappend A $Ar
        lappend b $br
    }
    set x [::math::linearalgebra::leastSquaresSVD $A $b]
} -result {1.0 1.0 0.1}

# Additional tests: procedures by Federico Ferri
#source ferri/ferri.test

set ::tcl_precision $old_precision

testsuiteCleanup