diff --git a/lib/la1.0/NIST.html b/lib/la1.0/NIST.html new file mode 100755 index 0000000000000000000000000000000000000000..11c744ca68e9868a3efcd04afa286d21959044e3 --- /dev/null +++ b/lib/la1.0/NIST.html @@ -0,0 +1,245 @@ + + +NIST PCA Example + +

+

NIST PCA Example

+
+# A Principal Components Analysis example
+#
+# walk through the NIST/Sematech PCA example data section 6.5
+# http://www.itl.nist.gov/div898/handbook/
+#
+# This file is part of the Hume La Package
+# (C)Copyright 2001, Hume Integration Software
+#
+# These calculations use the Hume Linear Algebra package, La
+#
+package require La
+catch { namespace import La::*}
+
+
+# Look at the bottom of this file to see the console output of this procedure.
+
+proc nist_pca {} {
+    puts {Matrix of observation data X[10,3]}
+    set X {2 10 3 7 4 3 4 1 8 6 3 5 8 6 1 8 5 7 7 2 9 5 3 3 9 5 8 7 4 5 8 2 2}
+    puts [show $X]
+    puts {Compute column norms}
+    mnorms_br X Xmeans Xsigmas	
+    puts "Means=[show $Xmeans]"
+    puts "Sigmas=[show $Xsigmas %.6f]"
+    puts {Normalize observation data creating Z[10,3]}
+    puts [show [set Z [mnormalize $X $Xmeans $Xsigmas]] %10.6f]\n
+    puts {Naive computation method that is sensitive to round-off error}
+    puts {ZZ is (Ztranspose)(Z) }
+    puts "ZZ=\n[show [set ZZ [mmult [transpose $Z] $Z]] %10.6f]"
+    puts {the correlation matrix, R = (Zt)(Z)/(n-1)}
+    puts "R=\n[show [set R [mscale $ZZ [expr 1/9.0]]] %8.4f]"
+    puts {Continuing the naive example use R for the eigensolution}
+    set V $R
+    mevsvd_br V evals
+    puts "eigenvalues=[show $evals %8.4f]"
+    puts "eigenvectors:\n[show $V %8.4f]\n"
+    puts "Good computation method- perform SVD of Z"
+    set U $Z
+    msvd_br U S V
+    puts "Singular values, S=[show $S]"
+    puts "square S and divide by (number of observations)-1 to get the eigenvalues"
+    set SS [mprod $S $S]
+    set evals [mscale $SS [expr 1.0/9.0]]
+    puts "eigenvalues=[show $evals %8.4f]"
+    puts "eigenvectors:\n[show $V %8.4f]"
+    puts "(the 3rd vector/column has inverted signs, thats ok the math still works)"
+    puts "\nverify V x transpose = I"
+    puts "As computed\n[show [mmult $V [transpose $V]] %20g]\n"
+    puts "Rounded off\n[show [mround [mmult $V [transpose $V]]]]"
+    #puts "Diagonal matrix of eigenvalues, L"
+    #set L [mdiag $evals]
+    #
+    
+    puts "\nsquare root of eigenvalues as a diagonal matrix"
+    proc sqrt x { return [expr sqrt($x)] }
+    set evals_sr [mat_unary_op $evals sqrt]
+    set Lhalf [mdiag $evals_sr]
+    puts [show $Lhalf %8.4f]
+
+    puts "Factor Structure Matrix S"
+    set S [mmult $V $Lhalf]
+    puts [show $S %8.4f]\n
+    set S2 [mrange $S 0 1]
+    puts "Use first two eigenvalues\n[show $S2 %8.4f]"
+    # Nash discusses better ways to compute this than multiplying SSt
+    # the NIST method is sensitive to round-off error
+    set S2S2 [mmult $S2 [transpose $S2]]
+    puts "SSt for first two components, communality is the diagonal"
+    puts [show $S2S2 %8.4f]\n
+
+
+    # define reciprocal square root function
+    proc recipsqrt x { return [expr 1.0/sqrt($x)] }
+    # use it for Lrhalf calculation
+    set Lrhalf [mdiag [mat_unary_op $evals recipsqrt]]
+
+    set B [mmult $V $Lrhalf]
+    puts "Factor score coefficients B Matrix:\n[show $B %12.4f]"
+
+    puts "NIST shows B with a mistake, -.18 for -1.2 (-1.18 with a typo)"
+
+    puts "\nWork with first normalized observation row"
+    set z1 {2 3 0 0.065621795889 0.316227766017 -0.7481995314}
+    puts [show $z1 %.4f]
+
+    puts "compute the \"factor scores\" from multiplying by B"
+    set t [mmult [transpose $z1] $B]
+    puts [show $t %.4f]
+    puts "NIST implies you might chart these"
+
+    #set T2 [dotprod $t $t]
+    #puts "Hotelling T2=[format %.4f $T2]"
+    
+    puts "Compute the scores using V, these are more familiar"
+    set t [mmult [transpose $z1] $V]
+    puts [show $t %.4f]
+    puts "Calculate T2 from the scores, sum ti*ti/evi"
+
+    set T2 [msum [transpose [mdiv [mprod $t $t] [transpose $evals]]]]
+    puts "Hotelling T2=[format %.4f $T2]"
+ 
+    puts "Fit z1 using first two principal components"
+    set p [mrange $V 0 1]
+    puts p=\n[show $p %10.4f]
+    set t [mmult [transpose $z1] $p]
+    set zhat [mmult $t [transpose $p]]
+    puts "  z1=[show $z1 %10.6f]"
+    puts "zhat=[show $zhat %10.6f]"
+    set diffs [msub [transpose $z1] $zhat]
+    puts "diff=[show $diffs %10.6f]"
+    puts "Q statistic - distance from the model measurement for the first observation"
+    set Q [dotprod $diffs $diffs]
+    puts Q=[format %.4f $Q]
+    puts "Some experts would compute a corrected Q that is larger since the 
+observation was used in building the model.  The Q calculation just used
+properly applies to the analysis of new observations."
+    set corr [expr 10.0/(10.0 - 2 -1)]  ;# N/(N - Ncomp - 1)
+    puts "Corrected Q=[format %.4f [expr $Q * $corr]] (factor=[format %.4f $corr])"
+    }
+
+return
+##########################################################################
+Console Printout
+##########################################################################
+Matrix of observation data X[10,3]
+7 4 3
+4 1 8
+6 3 5
+8 6 1
+8 5 7
+7 2 9
+5 3 3
+9 5 8
+7 4 5
+8 2 2
+Compute column norms
+Means=6.9 3.5 5.1
+Sigmas=1.523884 1.581139 2.806738
+Normalize observation data creating Z[10,3]
+  0.065622   0.316228  -0.748200
+ -1.903032  -1.581139   1.033228
+ -0.590596  -0.316228  -0.035629
+  0.721840   1.581139  -1.460771
+  0.721840   0.948683   0.676942
+  0.065622  -0.948683   1.389513
+ -1.246814  -0.316228  -0.748200
+  1.378058   0.948683   1.033228
+  0.065622   0.316228  -0.035629
+  0.721840  -0.948683  -1.104485
+
+Naive computation method that is sensitive to round-off error
+ZZ is (Ztranspose)(Z) 
+ZZ=
+  9.000000   6.017916  -0.911824
+  6.017916   9.000000  -2.591349
+ -0.911824  -2.591349   9.000000
+the correlation matrix, R = (Zt)(Z)/(n-1)
+R=
+  1.0000   0.6687  -0.1013
+  0.6687   1.0000  -0.2879
+ -0.1013  -0.2879   1.0000
+Continuing the naive example use R for the eigensolution
+eigenvalues=  1.7688   0.9271   0.3041
+eigenvectors:
+  0.6420   0.3847  -0.6632
+  0.6864   0.0971   0.7207
+ -0.3417   0.9179   0.2017
+
+Good computation method- perform SVD of Z
+Singular values, S=3.98985804571 2.88854344819 1.65449373616
+square S and divide by (number of observations)-1 to get the eigenvalues
+eigenvalues=  1.7688   0.9271   0.3041
+eigenvectors:
+  0.6420   0.3847   0.6632
+  0.6864   0.0971  -0.7207
+ -0.3417   0.9179  -0.2017
+(the 3rd vector/column has inverted signs, thats ok the math still works)
+
+verify V x transpose = I
+As computed
+                   1        -5.55112e-017         1.38778e-016
+       -5.55112e-017                    1         5.55112e-017
+        1.38778e-016         5.55112e-017                    1
+
+Rounded off
+1.0 0.0 0.0
+0.0 1.0 0.0
+0.0 0.0 1.0
+
+square root of eigenvalues as a diagonal matrix
+  1.3300   0.0000   0.0000
+  0.0000   0.9628   0.0000
+  0.0000   0.0000   0.5515
+Factor Structure Matrix S
+  0.8538   0.3704   0.3658
+  0.9128   0.0935  -0.3975
+ -0.4544   0.8838  -0.1112
+
+Use first two eigenvalues
+  0.8538   0.3704
+  0.9128   0.0935
+ -0.4544   0.8838
+SSt for first two components, communality is the diagonal
+  0.8662   0.8140  -0.0606
+  0.8140   0.8420  -0.3321
+ -0.0606  -0.3321   0.9876
+
+Factor score coefficients B Matrix:
+      0.4827       0.3995       1.2026
+      0.5161       0.1009      -1.3069
+     -0.2569       0.9533      -0.3657
+NIST shows B with a mistake, -.18 for -1.2 (-1.18 with a typo)
+
+Work with first normalized observation row
+0.0656 0.3162 -0.7482
+compute the "factor scores" from multiplying by B
+0.3871 -0.6552 -0.0608
+NIST implies you might chart these
+Compute the scores using V, these are more familiar
+0.5148 -0.6308 -0.0335
+Calculate T2 from the scores, sum ti*ti/evi
+Hotelling T2=0.5828
+Fit z1 using first two principal components
+p=
+    0.6420     0.3847
+    0.6864     0.0971
+   -0.3417     0.9179
+  z1=  0.065622   0.316228  -0.748200
+zhat=  0.087847   0.292075  -0.754958
+diff= -0.022225   0.024153   0.006758
+Q statistic - distance from the model measurement for the first observation
+Q=0.0011
+Some experts would compute a corrected Q that is larger since the 
+observation was used in building the model.  The Q calculation just used
+properly applies to the analysis of new observations.
+Corrected Q=0.0016 (factor=1.4286)
+
+ \ No newline at end of file diff --git a/lib/la1.0/NIST.tcl b/lib/la1.0/NIST.tcl new file mode 100755 index 0000000000000000000000000000000000000000..03439d5c60c4f03cc9fd79e07fcc330c9c60cee5 --- /dev/null +++ b/lib/la1.0/NIST.tcl @@ -0,0 +1,236 @@ +# A Principal Components Analysis example +# +# walk through the NIST/Sematech PCA example data section 6.5 +# http://www.itl.nist.gov/div898/handbook/ +# +# This file is part of the Hume La Package +# (C)Copyright 2001, Hume Integration Services +# +# These calculations use the Hume Linear Algebra package, La +# +package require La +catch { namespace import La::*} + + +# Look at the bottom of this file to see the console output of this procedure. + +proc nist_pca {} { + puts {Matrix of observation data X[10,3]} + set X {2 10 3 7 4 3 4 1 8 6 3 5 8 6 1 8 5 7 7 2 9 5 3 3 9 5 8 7 4 5 8 2 2} + puts [show $X] + puts {Compute column norms} + mnorms_br X Xmeans Xsigmas + puts "Means=[show $Xmeans]" + puts "Sigmas=[show $Xsigmas %.6f]" + puts {Normalize observation data creating Z[10,3]} + puts [show [set Z [mnormalize $X $Xmeans $Xsigmas]] %10.6f]\n + puts {Naive computation method that is sensitive to round-off error} + puts {ZZ is (Ztranspose)(Z) } + puts "ZZ=\n[show [set ZZ [mmult [transpose $Z] $Z]] %10.6f]" + puts {the correlation matrix, R = (Zt)(Z)/(n-1)} + puts "R=\n[show [set R [mscale $ZZ [expr 1/9.0]]] %8.4f]" + puts {Continuing the naive example use R for the eigensolution} + set V $R + mevsvd_br V evals + puts "eigenvalues=[show $evals %8.4f]" + puts "eigenvectors:\n[show $V %8.4f]\n" + puts "Good computation method- perform SVD of Z" + set U $Z + msvd_br U S V + puts "Singular values, S=[show $S]" + puts "square S and divide by (number of observations)-1 to get the eigenvalues" + set SS [mprod $S $S] + set evals [mscale $SS [expr 1.0/9.0]] + puts "eigenvalues=[show $evals %8.4f]" + puts "eigenvectors:\n[show $V %8.4f]" + puts "(the 3rd vector/column has inverted signs, thats ok the math still works)" + puts "\nverify V x transpose = I" + puts "As computed\n[show [mmult $V [transpose $V]] %20g]\n" + puts "Rounded off\n[show [mround [mmult $V [transpose $V]]]]" + #puts "Diagonal matrix of eigenvalues, L" + #set L [mdiag $evals] + # + + puts "\nsquare root of eigenvalues as a diagonal matrix" + proc sqrt x { return [expr sqrt($x)] } + set evals_sr [mat_unary_op $evals sqrt] + set Lhalf [mdiag $evals_sr] + puts [show $Lhalf %8.4f] + + puts "Factor Structure Matrix S" + set S [mmult $V $Lhalf] + puts [show $S %8.4f]\n + set S2 [mrange $S 0 1] + puts "Use first two eigenvalues\n[show $S2 %8.4f]" + # Nash discusses better ways to compute this than multiplying SSt + # the NIST method is sensitive to round-off error + set S2S2 [mmult $S2 [transpose $S2]] + puts "SSt for first two components, communality is the diagonal" + puts [show $S2S2 %8.4f]\n + + + # define reciprocal square root function + proc recipsqrt x { return [expr 1.0/sqrt($x)] } + # use it for Lrhalf calculation + set Lrhalf [mdiag [mat_unary_op $evals recipsqrt]] + + set B [mmult $V $Lrhalf] + puts "Factor score coefficients B Matrix:\n[show $B %12.4f]" + + puts "NIST shows B with a mistake, -.18 for -1.2 (-1.18 with a typo)" + + puts "\nWork with first normalized observation row" + set z1 {2 3 0 0.065621795889 0.316227766017 -0.7481995314} + puts [show $z1 %.4f] + + puts "compute the \"factor scores\" from multiplying by B" + set t [mmult [transpose $z1] $B] + puts [show $t %.4f] + puts "NIST implies you might chart these" + + #set T2 [dotprod $t $t] + #puts "Hotelling T2=[format %.4f $T2]" + + puts "Compute the scores using V, these are more familiar" + set t [mmult [transpose $z1] $V] + puts [show $t %.4f] + puts "Calculate T2 from the scores, sum ti*ti/evi" + + set T2 [msum [transpose [mdiv [mprod $t $t] [transpose $evals]]]] + puts "Hotelling T2=[format %.4f $T2]" + + puts "Fit z1 using first two principal components" + set p [mrange $V 0 1] + puts p=\n[show $p %10.4f] + set t [mmult [transpose $z1] $p] + set zhat [mmult $t [transpose $p]] + puts " z1=[show $z1 %10.6f]" + puts "zhat=[show $zhat %10.6f]" + set diffs [msub [transpose $z1] $zhat] + puts "diff=[show $diffs %10.6f]" + puts "Q statistic - distance from the model measurement for the first observation" + set Q [dotprod $diffs $diffs] + puts Q=[format %.4f $Q] + puts "Some experts would compute a corrected Q that is larger since the +observation was used in building the model. The Q calculation just used +properly applies to the analysis of new observations." + set corr [expr 10.0/(10.0 - 2 -1)] ;# N/(N - Ncomp - 1) + puts "Corrected Q=[format %.4f [expr $Q * $corr]] (factor=[format %.4f $corr])" + } + +return +########################################################################## +Console Printout +########################################################################## +Matrix of observation data X[10,3] +7 4 3 +4 1 8 +6 3 5 +8 6 1 +8 5 7 +7 2 9 +5 3 3 +9 5 8 +7 4 5 +8 2 2 +Compute column norms +Means=6.9 3.5 5.1 +Sigmas=1.523884 1.581139 2.806738 +Normalize observation data creating Z[10,3] + 0.065622 0.316228 -0.748200 + -1.903032 -1.581139 1.033228 + -0.590596 -0.316228 -0.035629 + 0.721840 1.581139 -1.460771 + 0.721840 0.948683 0.676942 + 0.065622 -0.948683 1.389513 + -1.246814 -0.316228 -0.748200 + 1.378058 0.948683 1.033228 + 0.065622 0.316228 -0.035629 + 0.721840 -0.948683 -1.104485 + +Naive computation method that is sensitive to round-off error +ZZ is (Ztranspose)(Z) +ZZ= + 9.000000 6.017916 -0.911824 + 6.017916 9.000000 -2.591349 + -0.911824 -2.591349 9.000000 +the correlation matrix, R = (Zt)(Z)/(n-1) +R= + 1.0000 0.6687 -0.1013 + 0.6687 1.0000 -0.2879 + -0.1013 -0.2879 1.0000 +Continuing the naive example use R for the eigensolution +eigenvalues= 1.7688 0.9271 0.3041 +eigenvectors: + 0.6420 0.3847 -0.6632 + 0.6864 0.0971 0.7207 + -0.3417 0.9179 0.2017 + +Good computation method- perform SVD of Z +Singular values, S=3.98985804571 2.88854344819 1.65449373616 +square S and divide by (number of observations)-1 to get the eigenvalues +eigenvalues= 1.7688 0.9271 0.3041 +eigenvectors: + 0.6420 0.3847 0.6632 + 0.6864 0.0971 -0.7207 + -0.3417 0.9179 -0.2017 +(the 3rd vector/column has inverted signs, thats ok the math still works) + +verify V x transpose = I +As computed + 1 -5.55112e-017 1.38778e-016 + -5.55112e-017 1 5.55112e-017 + 1.38778e-016 5.55112e-017 1 + +Rounded off +1.0 0.0 0.0 +0.0 1.0 0.0 +0.0 0.0 1.0 + +square root of eigenvalues as a diagonal matrix + 1.3300 0.0000 0.0000 + 0.0000 0.9628 0.0000 + 0.0000 0.0000 0.5515 +Factor Structure Matrix S + 0.8538 0.3704 0.3658 + 0.9128 0.0935 -0.3975 + -0.4544 0.8838 -0.1112 + +Use first two eigenvalues + 0.8538 0.3704 + 0.9128 0.0935 + -0.4544 0.8838 +SSt for first two components, communality is the diagonal + 0.8662 0.8140 -0.0606 + 0.8140 0.8420 -0.3321 + -0.0606 -0.3321 0.9876 + +Factor score coefficients B Matrix: + 0.4827 0.3995 1.2026 + 0.5161 0.1009 -1.3069 + -0.2569 0.9533 -0.3657 +NIST shows B with a mistake, -.18 for -1.2 (-1.18 with a typo) + +Work with first normalized observation row +0.0656 0.3162 -0.7482 +compute the "factor scores" from multiplying by B +0.3871 -0.6552 -0.0608 +NIST implies you might chart these +Compute the scores using V, these are more familiar +0.5148 -0.6308 -0.0335 +Calculate T2 from the scores, sum ti*ti/evi +Hotelling T2=0.5828 +Fit z1 using first two principal components +p= + 0.6420 0.3847 + 0.6864 0.0971 + -0.3417 0.9179 + z1= 0.065622 0.316228 -0.748200 +zhat= 0.087847 0.292075 -0.754958 +diff= -0.022225 0.024153 0.006758 +Q statistic - distance from the model measurement for the first observation +Q=0.0011 +Some experts would compute a corrected Q that is larger since the +observation was used in building the model. The Q calculation just used +properly applies to the analysis of new observations. +Corrected Q=0.0016 (factor=1.4286) \ No newline at end of file diff --git a/lib/la1.0/d2utmpYoxGTS b/lib/la1.0/d2utmpYoxGTS new file mode 100755 index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 diff --git a/lib/la1.0/defs b/lib/la1.0/defs new file mode 100755 index 0000000000000000000000000000000000000000..585baff24ca3898a7a0e55b44e0a60e239f75153 --- /dev/null +++ b/lib/la1.0/defs @@ -0,0 +1,92 @@ +# This file contains support code for the Tcl test suite. It is +# normally sourced by the individual files in the test suite before +# they run their tests. This improved approach to testing was designed +# and initially implemented by Mary Ann May-Pumphrey of Sun Microsystems. +# This file is NOT Copyright by Hume Integration +# +set VERBOSE 0 +set TESTS {} +set auto_noexec 1 + + +# If tests are being run as root, issue a warning message and set a +# variable to prevent some tests from running at all. + +set user {} +catch {set user [exec whoami]} + + +# Some of the tests don't work on some system configurations due to +# configuration quirks, not due to Tcl problems; in order to prevent +# false alarms, these tests are only run in the master source directory +# at Berkeley. The presence of a file "Berkeley" in this directory is +# used to indicate that these tests should be run. + +set atBerkeley [file exists Berkeley] + +proc print_verbose {test_name test_description contents_of_test code answer} { + puts stdout "\n" + puts stdout "==== $test_name $test_description" + puts stdout "==== Contents of test case:" + puts stdout "$contents_of_test" + if {$code != 0} { + if {$code == 1} { + puts stdout "==== Test generated error:" + puts stdout $answer + } elseif {$code == 2} { + puts stdout "==== Test generated return exception; result was:" + puts stdout $answer + } elseif {$code == 3} { + puts stdout "==== Test generated break exception" + } elseif {$code == 4} { + puts stdout "==== Test generated continue exception" + } else { + puts stdout "==== Test generated exception $code; message was:" + puts stdout $answer + } + } else { + puts stdout "==== Result was:" + puts stdout "$answer" + } +} + +proc test {test_name test_description contents_of_test passing_results} { + global VERBOSE + global TESTS + if {[string compare $TESTS ""] != 0} then { + set ok 0 + foreach test $TESTS { + if [string match $test $test_name] then { + set ok 1 + break + } + } + if !$ok then return + } + set code [catch {uplevel $contents_of_test} answer] + if {$code != 0} { + print_verbose $test_name $test_description $contents_of_test \ + $code $answer + } elseif {[string compare $answer $passing_results] == 0} then { + if $VERBOSE then { + print_verbose $test_name $test_description $contents_of_test \ + $code $answer + puts stdout "++++ $test_name PASSED" + } + } else { + print_verbose $test_name $test_description $contents_of_test $code \ + $answer + puts stdout "---- Result should have been:" + puts stdout "$passing_results" + puts stdout "---- $test_name FAILED" + } +} + +proc dotests {file args} { + global TESTS + set savedTests $TESTS + set TESTS $args + source $file + set TESTS $savedTests +} + diff --git a/lib/la1.0/index.html b/lib/la1.0/index.html new file mode 100755 index 0000000000000000000000000000000000000000..320081053e2824286ad61ab91a8742479895a71a --- /dev/null +++ b/lib/la1.0/index.html @@ -0,0 +1,10 @@ + + +La Package Documentation + + + + + + + diff --git a/lib/la1.0/la.html b/lib/la1.0/la.html new file mode 100755 index 0000000000000000000000000000000000000000..e18e0056a101ae85960978c85dcd365d35e625e0 --- /dev/null +++ b/lib/la1.0/la.html @@ -0,0 +1,1917 @@ + + + + + + + + +The Hume Linear Algebra Tcl Package + + + + + + +
+ +

The Hume Linear Algebra +Tcl Package, La

+ +

 

+ +

Feature Summary

+ +

 

+ +

The +package consists of Tcl procedures for the manipulation of vectors and +matrices.  The functionality spans +scaling, normalization, concatenation by rows or columns, subsetting by rows or +columns, formatted printing, transpose, dot product, matrix multiplication, +solution of linear equation sets, matrix inversion, eigenvalue/eigenvector +solutions,  singular value +decomposition, and solution of linear least squares.  The singular value decomposition can be used to perform the +principle components analysis of multivariate statistical process control and +avoid the possibly ill-conditioned multiplication of the XXt  matrix of  +observations.

+ +

 

+ +

The user can +mix vectors and arrays in linear algebra operations.  The logic does reasonable conversion of types.  Sophisticated operations such as evaluating +a custom procedure against each element of a matrix are easily possible.  Data is represented as ordinary Tcl list +variables so that the usual commands of Tcl are useable, and efficient access +to the data elements using compiled extensions is possible.

+ +

 

+ +

This +document and the package are ©Copyright 2001, Hume Integration Software.  You may use the package without licensing +fees according to the License Terms.

+ +

User +Guide

+ +

 

+ +

Obtaining +the Package

+ +

 

+ +

The +package may be obtained from the webpage http://www.hume.com/la.  +The package uses the "string is" Tcl command which was new +with Tcl 8.1.  So unless you are +comfortable making minor changes to the Tcl source code, you should only use +the package with Tcl version 8.1 and newer.  +The version 1.0 package only consists of text files, there are no binary +files.  The only difference between the +Windows archive and the POSIX archive is the difference in line feed / carriage +return characters at the end of each line of text.  The package is distributed as zip archives.  You are expected to use your own software to +unpack the archive files. 

+ +

Installation

+ +

 

+ +

The +package is installed by extracting from the zip archive file, the package +directory and files.  The package +directory should be added into the base directory of the Tcl runtime libraries.  So if your Tcl runtime library directory is +/usr/local/lib/tcl8.4, copy or extract the la1.0 directory and files to +/usr/local/lib so that you have the directory /usr/local/lib/la1.0 and the +package files are found as /usr/local/lib/la1.0/*.* .  The Tcl code should be compatible with Tcl version 8.1 and later +for any architecture, but we have only tested it with Tcl 8.3 and Tcl 8.4.   You can run the package regression tests by +navigating to the package directory and sourcing the file “la.test”.

+ +

 

+ +

Runtime Usage

+ +

 

+ +

The software +is found and made useable by your Tcl/Tk interpreter after you execute “package +require La”.  You typically import the +package procedure names into your global namespace to save yourself from having +to qualify the procedure names using the package name.  So typical initialization code looks like:

+ +

 

+ +

% package require La

+ +

1.1

+ +

% namespace import La::*

+ +

 

+ +

Operand Formats

+ +

The +package uses ordinary Tcl variables to represent scalars, vectors, and +arrays.  This means that the ordinary +Tcl commands are useable.  For example, +the Tcl set command can be used to copy a matrix:

+ +

 

+ +

set A [mident 5]  ;# A is assigned a 5x5 identity matrix

+ +

set B $A          ;# B is assigned the value of A

+ +

 

+ +

Tcl list +formats are used for better performance than arrays and for better +compatibility with C/C++ code.   +Dimension information is added in a simple way at the front of the list, +for vectors, matrices, and higher dimensional variables.

+ +

 

+ +

Scalar Format 

+ +

If the +list length of a variable is 1 the operand is a scalar, eg., +"6.02e23".  The Tcl command +llength can be used to determine the list length.

+ +

Vector Format

+ +

 

+ +

A vector of length N has the +representation as a Tcl list of:

+ +

 

+ +

2 N 0 v[0] v[1] v[2] ... v[N-1]

+ +

 

+ +

The first +element in the list, the value 2, signifies that there are two additional +elements of size information; the number of rows, and the number of +columns.  When the number of columns is +0, the operand is defined to be a one dimensional vector.  So the Tcl list sequence of {2 4 0 1 2 3 4} +represents the vector {1 2 3 4}.  A vector of length N has the same number of +elements as an Nx1 or 1xN matrix so they can be efficiently converted in place, +and indexing operations are simplified.

+ +

 

+ +

The index into the +underlying Tcl list for vector v[i] is

+ +

 

+ +

set index [expr {3 + $i}]

+ +

 

+ +

A vector +of length N is promoted to an N rows by 1 column matrix if used in an operation +where a matrix argument is expected.

+ +

 

+ +

The +transpose of a vector of length N is a 1 row by N columns matrix.

+ +

Matrix Format

+ +

 

+ +

The Tcl +list data used to represent two dimensional matrix a[R,C], has the format:

+ +

 

+ +

2 R C a[0,0] a[0,1] a[0,2] ... +a[0,C-1] \

+ +

a[1,0] a[1,1] ... a[1,C-1] \

+ +

... a[R-1,C-1] 

+ +

 

+ +

where R is +the number of rows, and C is the number of columns, and the backslash character +has been used to indicate that the multiple lines of text constitute a single +list sequence.  The format itself is +independent of whether the user thinks that indexing starts from 0 or 1.  In Tcl, indexing of lists always starts with +0 so we have chosen to consistently start indexing with 0.

+ +

 

+ +

For a +valid matrix, both the number of rows and the number of columns are positive +integers.  The data elements following +the {2 R C} dimension information can also be thought of as first row, second +row, …, last row.  The index into the +underlying Tcl list for a[i,j] is:

+ +

 

+ +

set index [expr {3 + $i * $C + +$j}]

+ +

 

+ +

Higher Dimensions

+ +

 

+ +

If the +first element in a Tcl list is > 2, the package assumes the list represents +a higher dimensional operand.  Logic for +higher dimension operands is not currently part of the package. 

+ +

 

+ +

The +candidate structure for 3-D data is:

+ +

 

+ +

3 P R C a[0,0,0] a[0,0,1] +a[0,0,2] ... a[0,0,C-1] ... a[P-1,R-1,C-1]

+ +

 

+ +

where +P=planes, R=rows, C=columns.  An +intuitive view is that the data is multiple 2-D planes of rows and columns +listed in order from the 0 plane to the P-1 plane. 

+ +

 

+ +

The index +into the underlying Tcl list for a[i,j,k] is

+ +

 

+ +

set index [expr 4 + $i*$R*$c + +$j*$c + $k]

+ +

 

+ +

Operand Examples

+ +

 

+ +

set pi 3.1415            ;# a scalar

+ +

set v {2 1 0 3.1415}    +;# v[1] with value 3.1415

+ +

# vectors should use the first dimension of 2

+ +

set v {2 0 1 3.1415}    +;# error - vector should use first dimension only

+ +

set v {2 3 0 1 2 3 4}   +;# error - vector has an additional element

+ +

set m {2 2 3 1 2 3 4 5 6}  ;# m[2,3]

+ +

show $m

+ +

1 2 3

+ +

4 5 6

+ +

show [transpose $m]

+ +

1 4

+ +

2 5

+ +

3 6

+ +

 

+ +

Argument Passing By Value and By Reference

+ +

There are +typically two procedures defined for an algorithm.  A plainly named procedure such as "transpose", and +another procedure with the suffix _br such as +"transpose_br".  The plainly +named procedures expect their data arguments to be passed by value, which is +the usual argument passing convention with Tcl programming.  The plain calls are designed for ease of +interactive use, and in general have been coded to perform more conversion of +arguments, more error checking, and to trade off efficiency for +convenience.  The _br procedures +are intended for efficent use, with the _br indicating that data +arguments are passed "By Reference" to avoid copying the data.  In Tcl, to pass by reference means that the +caller has the data in a named variable, and the caller passes the name of the +variable instead of a copy of the data.  +You can see that passing by reference is more efficient for larger +vectors and arrays.  The _br +procedures in general assume that the data arguments have the correct +structure, so the caller may need to use the promote, demote, and +transpose calls to prepare arguments for a _br call.

+ +

 

+ +

A +Note on Performance

+ +

 

+ +

Number +crunching in Tcl?  The conventional +wisdom is that computationally intensive tasks such as numerical analysis are +high on the list of things that you should not do in a scripting language.   It is true that you can obtain better +performance using a traditional compiled language such as Fortran, or C +code.   But nowadays, the performance +gain that you will obtain by switching languages is similar to the performance +gain you will see between a new computer and a computer from two years +ago.   In many situations, the actual +time spent processing numbers is a small portion of the overall project.  The time it takes the researcher to explore +his data, decide on the desired approach, and develop the software for his +desired analysis is the only performance time that really matters.   Here, a high-level development environment +such as Tcl or MATLAB wins over C++ or Java hands down.   Actual performance numbers for the Tcl +package are nothing to be ashamed about.   +Using a 900Mhz Pentium III notebook computer with Tcl 8.4.1, the following times +were obtained for inverting an NxN matrix using Gauss elimination with partial +pivoting.   This is something of a worst +case since most numerical algorithms avoid the direct inversion of a full +matrix.  Often a major performance +improvement is made by choosing a better algorithm.

+ +

 

+ +

 

+ + + + + + + + + + + + + + + + + + + + + + + + + +
+

Inversion + of an NxN matrix.

+

 

+

Development + time for the following timed benchmark: 15 seconds.  It worked properly on the first attempt.  The flush and update commands are not + needed except to show progress before completion.  The benchmark:

+

 

+

% source la.tcl

+

% namespace import La::*

+

% foreach n {10 20 50 100} {

+

   puts [time + {msolve [mdingdong $n] [mident $n]}]

+

   flush stdout

+

   update

+

 }

+
+

Matrix + Dimension (N)

+
+

Inversion + Time (seconds)

+
+

10

+
+

0.012

+
+

20

+
+

0.085

+
+

50

+
+

1.251

+
+

100

+
+

9.965

+
+ +

 

+ +

Here are +some tips for efficient coding.  In Tcl +the expr command is used for the evaluation of +numeric expressions.  You should almost +always surround the arguments to expr with braces to prevent the Tcl +interpreter from substituting the arguments.  +The expr command itself is able to substitute variable references +directly and avoid re-interpreting their string representation.  For example:

+ +

 

+ +

set index [expr {3 + $i * $C + +$j}]   ;# good example

+ +

set index [expr 3 + $i*$C + +$j]       ;# less efficient example

+ +

 

+ +

The +Windows NT/Windows 2000 Tk console does not perform very well when large text +strings are displayed.  You may find it +preferable to print results by separate rows for large output sets, or to write +large output results to a file.  Another +alternative to the Windows Tk console, is to use tclsh83.exe from a command +window.

+ +

 

+ +

Tcl +maintains an internal binary representation for variable values, and computes a +string representation for a variable value only when needed.   If your logic is inconsistent and sometimes +treats a variable as a string, and sometimes treats it as a list of numbers, +you can cause inefficient “shimmering” of the internal representation.  For best efficiency use list oriented +commands to manipulate your vector and matrix variables, and don’t mix string +oriented commands.  For example:

+ +

 

+ +

      lappend m +$newvalue      ;# good – lappend is a list +oriented procedure

+ +

      append m +" $newvalue"    ;# bad, this +causes conversion to a string

+ +

 

+ +

The last +tip is that it is very easy to use the built-in time +command to analyze actual performance of your code. 

+ +

 

+ +

A PCA +Example

+ +

 

+ +

The +package contains a worked Principle Components Analysis problem based on +Section 6.5 of the SEMATECH/NIST Statistics Handbook, http://www.itl.nist.gov/div898/handbook.

+ +

 

+ +

See the +file NIST.tcl.

+ +

Reference

+ +

 

+ +

Alphabetical Procedure +Table

+ +

 

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+

Procedure + Name(s)

+
+

Description

+
+

demote, demote_br

+
+

Demote + an Nx1 or 1xN matrix to a vector[N].  + Demote a vector[1] to a scalar.  + Call twice to demote a 1x1 matrix to scalar.

+
+

dim, dim_br

+
+

Return + the dimension of the argument,  and to + some degree verify a proper format.

+
+

dotprod, dotprod_br

+
+

Dot Product + = sum over i, Ai * Bi

+

Can work + columns or rows in matrices because indexing increments are optional + arguments.

+
+

join_cols, join_cols_br

+
+

Combine + vector or matrices as added columns.

+
+

join_rows, join_rows_br

+
+

Combine + vector or matrices as added rows.

+
+

lassign, lassign_br

+
+

Replace + a single value in a Tcl list.

+
+

madd

+
+

Computes + a new matrix or vector from the addition of corresponding elements in two + operands.

+
+

madjust

+
+

Apply a + scale factor and offset to the operand’s elements.

+
+

mat_binary_op, mat_binary_op_br

+
+

Perform + binary operations such as addition on corresponding elements of vectors or + matrices.

+
+

mat_unary_op, mat_unary_op_br

+
+

Perform + unary operations on elements like scaling.

+
+

mathprec

+
+

Returns + the smallest number epsilon, such that 1+epsilon > 1.

+

Also + reports on the machine radix, digits, and rounding behavior

+
+

mcols, mcols_br

+
+

Return + the number of columns in a matrix or vector.  + You can use mrows and mcols to keep your software isolated from the + details of the data representation.

+
+

mdiag

+
+

Creates + a diagonal matrix from a vector, an Nx1 matrix, or a 1XN matrix.

+
+

mdingdong

+
+

Create + the Ding Dong test matrix, a Cauchy matrix that is represented inexactly in + the machine, but very stable for inversion by elimination methods.

+
+

mdiv

+
+

Computes + a new matrix or vector from the division of corresponding elements in two + operands.

+
+

mevsvd, mevsvd_br

+
+

Solve + for the eigenvectors and eigenvalues of a real symmetric matrix by singular + value decomposition.

+
+

mhilbert

+
+

Create + the Hilbert test matrix which is notorious for being  ill conditioned for eigenvector/eigenvalue + solutions.

+
+

mident

+
+

Create + an identity matrix of order N.

+
+

mlssvd

+
+

Linear + least squares solution for over-determined linear equations using singular + value decomposition.

+
+

mmult, mmult_br

+
+

Ordinary + matrix multiplication.

+
+

mnorms, mnorms_br

+
+

Compute + the means and standard deviations of each column.

+
+

mnormalize, mnormalize_br

+
+

Normalize + each column by subtracting the corresponding mean and then dividing by the + corresponding standard deviation.

+
+

moffset

+
+

Add a + constant to the elements of an operand.

+
+

mprod

+
+

Computes + a new matrix or vector from the product of corresponding elements in two + operands.

+
+

mrange, mrange_br

+
+

Return a + subset of selected columns, selected rows as a new matrix.  Also can be used to reverse the ordering + when the start index > end index.

+
+

mround, mround_br

+
+

Round + off elements in a matrix if they are close to integers.

+
+

mrows, mrows_br

+
+

Return + the number of rows in a matrix or vector.  + You can use mrows and mcols to keep your software isolated from the + details of the data representation.

+
+

mscale

+
+

Multiply + each element in an operand by a factor.

+
+

msolve, msolve_br

+
+

Solve + the matrix problem Ax = p for x, where p may be multiple columns.  When p is the identity matrix, the + solution x, is the inverse of A.  Uses + Gauss elimination with partial pivoting.

+
+

msub

+
+

Computes + a new matrix or vector from the subtraction of corresponding elements in two + operands.

+
+

msum, msum_br

+
+

Compute + the sums of each column, returning a vector or scalar result.  Call twice to get the total sum of columns + and rows (set total [msum [msum $a]]).

+
+

msvd

+
+

Perform + the Singular Value Decomposition of a matrix.

+
+

promote, promote_br

+
+

Promote + a scalar or vector to an array.  + Vector[N] is promoted to an Nx1 array.

+
+

show, show_br

+
+

Return a + formatted string representation for an operand.  Options allow for specify the format of numbers, and the + strings used to separate column and row elements.

+
+

transpose, transpose_br

+
+

Performs + the matrix transpose, exchanging [i,j] with [j,i].  A vector is promoted to a 1xN array by transpose.

+
+

vdiag, vdiag_br

+
+

Create a + vector from the diagonal elements of a matrix.

+
+

vtrim, vtrim_br

+
+

For a + vector or matrix operand, just return the actual data elements by trimming + away the dimension and size data in the front

+
+ +

 

+ +

 

+ +

demote x

+ +

demote_br name_x_in +{name_out {}}

+ +

 

+ +

Demote an Nx1 or 1xN matrix to a vector[n].  Demote a vector[1] to a scalar.  Call twice to demote a 1x1 matrix to +scalar.  For the _br procedure, the +default output destination is to overwrite the input.

+ +

 

+ +

dim x

+ +

dim_br name_x

+ +

 

+ +

Return +the dimension of the argument,  and to +some degree verify a proper format.  Returns +0 for a scalar, 1 for a vector, 2 for a matrix, and an empty string, {}, for an +empty input value.  Most improper +formats will result in an error.

+ +

dotprod a b {N {}} {a0index 3} {b0index 3} {ainc 1} +{binc 1}

+ +

dotprod_br a_name b_name {N {}} {a0index 3} {b0index 3} {ainc 1} {binc 1}

+ +

Perform the dot product of two operands, a +and b, returning a scalar result.  +The default arguments will correctly process vector and conformable 1xN +or Nx1 matrices.  The argument N +is the number of element pairs to multiply when computing the sum of pair +products.  The a0index and b0index +optional arguments are the offset of the first element in the a and b +operands, respectively.  The ainc +and binc optional arguments are the index increment values to be added +to the indexes into a and b to obtain subsequent elements.

+ +

 

+ +

join_cols a b

+ +

join_cols_br {a_in b_in +{c_out {}}

+ +

Combine vector or matrices as added columns, +returning a matrix result.  The +arguments must have the same number of rows.  +The default output for the join_cols_br call is to overwrite the +a input.

+ +

 

+ +

join_rows a b

+ +

join_rows_br {a_in b_in +{c_out {}}

+ +

Combine vector or matrices as added rows, returning +a matrix result.  The arguments must +have the same number of columns.  The +default output for the join_rows_br call is to overwrite the a +input.

+ +

lassign list index value

+ +

lassign_br listname index +value

+ +

Replace a single element in a Tcl list.  There is conditional logic in the package so that the lassign_br +could be replaced by a C code version that would be able to update the list +directly without copying the data the way lreplace +does.

+ +

madd a b

+ +

This procedure uses mat_binary_op to return a matrix or +vector result from the addition of corresponding elements in two operands.

+ +

madjust a scale offset

+ +

This procedure uses mat_unary_op to return a matrix or +vector where the elements have been multiplied by a common factor, scale, +and then the offset value has been added.

+ +

 

+ +

mat_binary_op a b {op +}

+ +

mat_binary_op_br a_in +b_in op {c_out {}}

+ +

This procedure is used to execute binary operations +such as addition on corresponding elements of vector or matrix operands.  The default output for the _br call is to +overwrite the a input.

+ +

mathprec {puts puts}

+ +

Returns the smallest number epsilon, such that +1+epsilon > 1.  Also reports on the +machine radix, digits, and rounding behavior, by using $puts as a command with +a single string argument, of the format “radix=2.0 digits=53 +epsilon=2.22044604925e-016 method=truncation”.  +The default value of the puts argument causes the string to be printed +at the console.

+ +

mcols m

+ +

mcols_br m_name

+ +

Return the number of columns in a matrix or +vector.  You can use mrows +and mcols to keep your software isolated from the details of the data +representation.

+ +

 

+ +

mdiag v

+ +

Creates a diagonal matrix from a vector, an Nx1 matrix, or a 1xN +matrix.

+ +

mdingdong N

+ +

Create the Ding Dong test matrix of size NxN.  The matrix is a Cauchy matrix that is represented +inexactly in the machine, but very stable for inversion by elimination +methods.  Created by Dr. F. N. Ris +(Nash, 1979).

+ +

mdiv a b

+ +

Uses mat_binary_op to compute a new matrix or vector from +the division of corresponding elements in two operands. 

+ +

mevsvd A {epsilon 2.3e-16}

+ +

mevsvd_br A_in_out +evals_out {epsilon 2.3e-16}

+ +

Solve for the eigenvectors and eigenvalues of a real symmetric +matrix by singular value decomposition.  +The eigenvectors of the solution are returned as the columns of A.  The epsilon argument is expected to be the +value returned by the mathprec procedure for your platform.

+ +

mhilbert N

+ +

Create the Hilbert test matrix which is notorious for being  ill conditioned for eigenvector/eigenvalue +solutions. (Nash, 1979)

+ +

mident N

+ +

Create an identity matrix of order N.

+ +

mlssvd A y {q 0.0} {puts puts} {epsilon +2.3e-16}

+ +

Solve the linear least squares solution for over-determined linear +equations using singular value decomposition.  +Solves the problem A[m,n]x[n] ~ y[m] for x[n] +where each row of A is a set of dependent variable values, x[n] +is the vector of independent variables, and y[m] is the set of dependent +values such as measured outcome values.  +The first column of A is usually all ones to compute a constant +term in the regression equation.  The value +q is specified such that singular values less than q are treated +as zero.  Typically a judgment is made +as to what variation is significant, and what is just noise.  The significance level varies by +application.  For example, a small value +of q might be appropriate to detect a new planet from orbital data, +whereas larger values would be used to regress socio-economic statistics.  The default value of the puts +argument causes the singular values to be printed at the console after the +matrix is factored.  The epsilon +argument is expected to be the result of the mathprec procedure for your +platform.

+ +

mmult A B

+ +

mmult_br name_A name_B +C_out

+ +

Ordinary matrix multiplication, A[p,q] x B[q,r] = C[p,r].  Vector arguments are  promoted to Nx1 arrays, so chances are if +you are using one as a left operand you probably intend to use the transpose of +it (1xN), which is easily done using the transpose +procedure.  The mmult procedure +returns the value of the matrix product.  +The mmult_br procedure writes the matrix product into the +variable whose name is C_out.  +The variable name C_out should specify a different variable than +either of the input variables specified by name_A or name_B.

+ +

mnorms a

+ +

mnorms_br name_a +means_out sigmas_out

+ +

Compute the means and standard deviations of each column of the +input matrix.  Vector results are +returned.  If the number of rows is less +than 2, the standard deviations cannot be calculated and an error is returned.

+ +

mnormalize a means sigmas

+ +

mnormalize_br name_a means_in +sigmas_in {c_out {}}

+ +

Normalize each column of a matrix by subtracting the corresponding +mean and then dividing by the corresponding standard deviation.  The default output matrix for mnormalize_br +is to overwrite the input matrix specified by name_a.

+ +

moffset a delta

+ +

Uses mat_unary_op to add a scalar constant to the elements +of an operand.  The modified operand is +the returned result.

+ +

mprod a b

+ +

Uses mat_binary_op to compute a new matrix or vector from the +multiplicative product of corresponding elements in two operands.

+ +

mrange m col_start col_last {row_start 0} +{row_last end}

+ +

mrange_br name_m_in c_out +col_start col_last {row_start 0} {row_last end}

+ +

Returns a subset of the selected columns and selected rows of a +matrix as a new matrix.  Column and row +indexing begin with 0.  The token end +may be used to specify the last row or column.  The default values for row selection will return all rows.  The rows and columns are read and copied +from the input matrix in the order of start to last.  The selection result includes rows and columns with index values +equal to the start indexes and proceeds to the last indexes, including all rows +and columns with indexes between the start and last values and equal to the start +and last values.  If the start index is +less than the last index for row or column selection, the result matrix is +created with the row or columns in reverse order.  This is a feature, not a bug.  +

+ +

% set m {2 2 3 1 2 3 4 5 6}

+ +

% show $m

+ +

1 2 3

+ +

4 5 6

+ +

% show [mrange $m 0 1]

+ +

1 2

+ +

4 5

+ +

% show [mrange $m end 0]

+ +

3 2 1

+ +

6 5 4

+ +

% show [mrange $m 1 1 0 0]

+ +

2

+ +

mround a {epsilon 1.0e-8}

+ +

Uses mat_unary_op to round-off to the nearest +integer the elements of a matrix or vector which are within epsilon of an +integer value.

+ +

% show $m %12.4g

+ +

           1   +1.388e-016   8.327e-017  -2.776e-017

+ +

  1.388e-016            +1  -5.551e-017   -1.11e-016

+ +

  8.327e-017  +-5.551e-017            1  -5.551e-017

+ +

 -2.776e-017   +-1.11e-016  -5.551e-017            1

+ +

% show [mround $m]

+ +

1.0 0.0 0.0 0.0

+ +

0.0 1.0 0.0 0.0

+ +

0.0 0.0 1.0 0.0

+ +

0.0 0.0 0.0 1.0

+ +

mrows +m

+ +

mrows_br m_name

+ +

Return the number +of rows in a matrix or vector.  You can +use mrows and mcols to keep +your software isolated from the details of the data representation.

+ +

mscale a scalefactor

+ +

Uses mat_unary_op to multiply each element in +a matrix or vector operand, a, by the scalar scalefactor.

+ +

msolve A p

+ +

msolve_br Ap_in {tolerance 2.3e-16}

+ +

Solves a system of linear +equations, Ax = p, for x, by brute force application of +Gauss elimination with partial pivoting.  +When p is the +identity matrix, the solution x, is the inverse of matrix A.The msolve_br procedure accepts as input the +columnwise concatenation of A and p, and overwrites the p +columns with the solution columns x.   +The msolve procedure accepts A and p as separate +arguments, and returns the solution data directly.  The epsilon +argument is expected to be the result of the mathprec procedure for your +platform.

+ +

msub a b

+ +

Uses mat_binary_op to compute a new matrix or vector from +the subtraction of corresponding elements in two operands.

+ +

msum a

+ +

msum_br name_a sums_out

+ +

Compute the sums of each column of a matrix or +vector, returning a vector or scalar result.  +Call twice to get the total sum of columns and rows (set total [msum +[msum $a]]).

+ +

msvd a

+ +

msvd_br name_a_in_U_out S_out V_out {epsilon 2.3e-16}

+ +

Perform the Singular Value Decomposition of a +matrix.

+ +

This factors matrix A into (U)(S)(Vtrans) where

+ +

    A[m,n] is the original +matrix

+ +

    U[m,n] has orthogonal columns (Ut)(U) =   (1(k)

+ +

       and multiplies to an identity matrix     ...

+ +

       supplemented with zeroes if needed        0(n-k))

+ +

    V[n,n] is orthogonal   +(V)(Vtran) = [mident $n]

+ +

       V contains the eigenvectors aka the principal components

+ +

    S is diagonal with the positive singular values of A

+ +

       Square S and divide by (m-1) to get the principal component

+ +

       eigenvalues.

+ +

   

+ +

    A[m,n]V[n,n] = B[m,n]  transforms +A to orthogonal columns, B

+ +

    B[m,n] = U[m,n]S[n,n] 

+ +

   

+ +

The msvd procedure outputs +formatted results to the console.  The msvd_br +procedure overwrites the input matrix a with the U result +matrix.  The epsilon argument is expected to be the result of the mathprec +procedure for your platform.

+ +

mat_unary_op a op

+ +

mat_unary_op_br name_a op {name_c_out {}}

+ +

Perform unary operations on operand elements like scaling.  The default output of the mat_unary_op_br +procedure is to overwrite the input specified by name_a.  The logic sweeps through the matrix or +vector, and evaluates the concatenation of $op with the operand +element.  For example, if the value of op +is “expr 0.5 *” the result would be to multiply each element by 0.5.  You can define your own procedures, and pass +the procedure name as the op argument.

+ +

promote x

+ +

promote_br name_x {name_out {}}

+ +

Promote a scalar or vector to an array.  Vector[N] is promoted to an Nx1 array.  Calling promote with a matrix argument is all right; the result +is the matrix unchanged.  The default +output of promote_br is to overwrite the input.

+ +

show x {format {}} {col_join { }} +{row_join \n}}

+ +

show_br name_in {name_out +{}} {format {}} {col_join { }} {row_join \n}

+ +

Return a formatted string representation for an operand.   Options allow for specify the format of +numbers, and the strings used to separate column and row elements.  The optional format argument is a +specification string which is applied to convert each element using the Tcl +format command.  +The format command is very is similar to the C code sprintf +command.  

+ +

 

+ +

set m {2 2 3 1 2 3 4 5 6}

+ +

% show $m %6.2f

+ +

  1.00   2.00   3.00

+ +

  4.00   5.00   6.00

+ +

% show $m {} , \;

+ +

1,2,3;4,5,6

+ +

transpose x

+ +

transpose_br name_x +{name_out {}}

+ +

Performs the matrix transpose, exchanging [i,j] with +[j,i].  A vector is promoted to a 1xN +array by transpose.  The default output +of the transpose_br procedure is to overwrite the input data named by name_x.

+ +

vdiag m

+ +

vdiag_br m_name_in v_out

+ +

Create a vector from the diagonal elements of a +matrix.

+ +

vtrim x

+ +

vtrim_br name_x {name_out +{}}

+ +

For a vector or matrix operand, just return the +actual data elements by trimming away the dimension and size data in the front +of the underlying Tcl list representation.  +The default output of the vtrim_br procedure is to overwrite the +input data named by name_x.

+ +

About the Author

+ +

This package has been developed by Edward C. Hume, III +PhD.  Dr. Hume has been interested in +numerical methods since the early 1980’s when his doctoral research at MIT +involved a comparison of Finite Element and Boundary Element methods for moving +boundary problems.  In recent years he +has been applying univariate and multivariate Statistical Process Control +techniques in his consulting work.  Dr. +Hume is the founder of Hume Integration Software, a software product and +consulting company, with a focus on Computer Integrated Manufacturing in the +Semiconductor and Electronics industries.

+ +

 

+ +

Hume Integration's flagship product is the Distributed Message +Hub (DMH) Application Development Package, a cohesive and synergistic set +of tools that extends the Tcl/Tk programming environment. The package includes +an in-memory SQL database that has subscription capability, comprehensive +support for equipment interfaces using SECS, serial, or network protocols, and +high-level facilities for interprocess communication. Applications can easily +share data and logic by exchanging Tcl and SQL messages which are efficiently +processed by the extended interpreter. This toolset is offered for Windows +2000/NT and major POSIX platforms including HP-UX, Linux, AIX, and SOLARIS. The +toolset is in 7x24 use in dozens of factories located in the United States, +Malaysia, Korea, Japan, China, Hong Kong, Singapore, Taiwan, Mexico, France, +and Scotland.  Feel free to visit our +website at http://www.hume.com.

+ +

 

+ +

References

+ +

 

+ +

The more +sophisticated algorithms in this package were adapted from Nash, 1979.

+ +

 

+ +

Compact +Numerical Methods for Computers: Linear Algebra and Function Minimisation by J. C. Nash, John Wiley & +Sons, New York, 1979.

+ +

License +Terms

+ +

 

+ +

The La +package software is being distributed under terms and conditions similar to +Tcl/Tk.  The author is providing the +package software to the Tcl community as a returned favor for the value of +packages received over the years.

+ +

 

+ +

The La package software is copyrighted by Hume Integration +Software. The following terms apply to all files associated with the software +unless explicitly disclaimed in individual files.

+ +

 

+ +

The authors hereby grant permission to use, copy, modify, +distribute, and license this software and its documentation for any purpose, +provided that existing copyright notices are retained in all copies and that +this notice is included verbatim in any distributions. No written agreement, +license, or royalty fee is required for any of the authorized uses. +Modifications to this software may be copyrighted by their authors and need not +follow the licensing terms described here, provided that the new terms are +clearly indicated on the first page of each file where they apply.

+ +

 

+ +

IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO +ANY PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES +ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES +THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGE.

+ +

 

+ +

THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY +WARRANTIES,INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN "AS +IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE +MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

+ +

 

+ +

GOVERNMENT USE: If you are acquiring this software on behalf +of the U.S. government, the Government shall have only "Restricted +Rights" in the software and related documentation as defined in the +Federal Acquisition Regulations (FARs) in Clause 52.227.19 (c) (2).  If you are acquiring the software on behalf +of the Department of Defense, the software shall be classified as +"Commercial Computer Software" and the Government shall have only +"Restricted Rights" as defined in Clause 252.227-7013 (c) (1) of +DFARs.  Notwithstanding the foregoing, +the authors grant the U.S. Government and others acting in its behalf +permission to use and distribute the software in accordance with the terms +specified in this license.

+ +

 

+ +

Document +Version

+ +

Date of +last revsion $Date: 2002/11/05 19:45:52 $.

+ +
+ + + + diff --git a/lib/la1.0/la.tcl b/lib/la1.0/la.tcl new file mode 100755 index 0000000000000000000000000000000000000000..f9b6e19e386e053bcb4ea52d8c5d71d06f0fd51d --- /dev/null +++ b/lib/la1.0/la.tcl @@ -0,0 +1,1549 @@ +# $Header: /usr/cvsroot/la/la.tcl,v 1.4 2002/11/05 19:45:52 hume Exp $ +# +# La == Linear Algebra and similar math +# +# Part of the La Package +# (C)Copyright 2001, Hume Integration Software +# +# $Log: la.tcl,v $ +# Revision 1.4 2002/11/05 19:45:52 hume +# Updated for Tcl 8.4, fixed mround procedure, found 8.4 lset slower than +# original lassignLa_br code. +# +# Revision 1.3 2001/11/14 13:27:35 hume +# Changed Services to Software. +# +# Revision 1.2 2001/07/12 23:59:51 hume +# Added K proc and boosted performance bigtime! +# +# Revision 1.1.1.1 2001/07/12 15:50:16 hume +# Initial check-in. +# +# +# + +package provide La 1.0.1 + + + +# +# Package scope and addressed requirements: +# + +# Procedures cover manipulation of vectors and matrices such as +# scaling, concatenation by rows or columns, subsetting by +# rows or columns, formatted printing, transpose, dot product, +# matrix multiplication, solution of linear equation sets, +# matrix inversion, and eigenvalue/eigenvector solutions. +# The user can mix vectors and arrays in linear algebra operations. +# The logic does reasonable conversion of types. +# Sophisticated operations such as evaluating a custom procedure +# against each element of a matrix are easily possible. +# +# There are typically two calls for a procedure. A plainly +# named procedure such as "transpose", and another procedure +# with the suffix "_br" such as "transpose_br". The plainly +# named procedures expect their data arguments to be passed +# by value, which is the usual argument passing convention with +# Tcl programming. The plain calls are designed for ease of +# interactive use, and in general have been coded to perform +# more conversion of arguments, more error checking, and to +# trade off efficiency for convenience. The "_br" procedures +# are intended for efficent use, with the "_br" indicating that +# data arguments are passed "By Reference" to avoid copying the +# data. In Tcl, to pass by reference means that the caller +# has the data in a named variable, and the caller passes the +# name of the variable instead of a copy of the data. You can +# see that passing by reference becomes important for larger +# vectors and arrays. The _br procedures in general assume +# that the data arguments have the correct structure, so the +# caller may need to use the "promote", "demote", and "transpose" +# calls to prepare arguments for a _br call. + + + +# +# La operand formats - scalars, vectors, and arrays +# +# Tcl list formats are used for better performance than arrays +# and better compatibility with C/C++ code. +# We add dimension information in a simple way at the front of the list. +# +# La operands are structured Tcl lists. +# If the list length is 1 the operand is a scalar, eg., "6.02e23". +# If the first element in the list is 2, the two following +# elements provide the number of rows and the number of columns. +# If the first element is 2: +# If the number of columns is 0, the operand is a one dimensional +# vector. So {2 4 0 1 2 3 4} is the vector {1 2 3 4}. Vectors +# have the same structure as matrices so they can be efficiently +# modified in place. +# So a vector of length N has the representation: +# {2 N 0 v0 v1 v2 ... vN-1} +# The index into the underlying Tcl list for v[i] is +# set index [expr 3 + $i] +# A vector of length N is promoted to an N rows by 1 column +# matrix if used in an operation where a matrix argument is +# expected. +# The transpose of a vector of length N is a 1 row by N columns +# matrix. +# If the number of rows and the number of columns are both positive, +# the operand is a matrix. The data elements follow the number of +# rows and columns. The data is concatenated as row0, row1, ... +# so a 2-D matrix, a, has the representation: +# 2 R C a[0,0] a[0,1] a[0,2] ... a[0,C-1] a[1,0] a[1,1] ... a[1,C-1] ... a[R-1,C-1] +# +# The index into the underlying Tcl list for a[i,j] is +# set index [expr 3 + $i * $C + $j] +# +# If the first element in the list is > 2, the package assumes the list +# represents a higher dimensional operand. Logic for higher dimension +# operands is not currently part of the package. +# The candidate structure for 3-D data is: +# 3 P R C a[0,0,0] a[0,0,1] a[0,0,2] ... a[0,0,C-1] ... a[P-1,R-1,C-1] +# P=planes R=rows, C=cols +# An intuitive view is that the data is multiple 2-D planes of rows +# and columns listed in order from the 0 plane to the P-1 plane. +# The index into the underlying Tcl list for a[i,j,k] is +# set index [expr 4 + $i*$R*$c + $j*$c + $k] +# +# +# +# La operands, explained by example: +# +# a one element scalar (llength == 1) +# eg., 3.1415 +# vectors should use the first dimension of 2 +# 2 1 0 3.1415 ;# a vector of length 1 +# 2 0 1 3.1415 ;# error - vector should use first dimension only +# 2 3 0 1 2 3 ;# vector {1 2 3} +# 2 3 0 1 2 3 4 ;# error - vector has additional element +# 2 0 0 ;# error - cannot distinguish 1-D or 2-D +# 2 0 0 ;# error - simpler to have scalars always llength 1 +# a 2D array of N rows, M columns +# 2 N M a[0,0] a[0,1] a[0,2] ... a[0,M] a[1,0] a[1,1] ... a[1,M] ... a[N,M] +# + +# getting started +# +# % source "la.tcl" ;# advanced users will use "package require La" +# % namespace import La::* ;# makes names like "La::show" useable as "show" +# set x {2 2 3 1 2 3 4 5 6} +# % show $x +# 1 2 3 +# 4 5 6 +# % show [transpose $x] +# 1 4 +# 2 5 +# 3 6 +# % set tx [transpose $x] +# 2 3 2 1 4 2 5 3 6 +# % show [join_cols $tx $tx] +# 1 4 1 4 +# 2 5 2 5 +# 3 6 3 6 +# % show [moffset $x 0.2] +# 1.2 2.2 3.2 +# 4.2 5.2 6.2 +# % set v {2 3 0 1 2 3} +# % show [mmult $x $v] +# 14.0 +# 32.0 + + + +############################### Package Source Code ##################### + +# Advanced Users: +# The package namespace name is defined here by default. +# You can specify a different name before using the package. +global La_namespace +if { ![info exists La_namespace] } { set La_namespace La } + +namespace eval $La_namespace { + + namespace export \ +dim dim_br \ +dotprod dotprod_br \ +join_cols join_cols_br \ +join_rows join_rows_br\ +lassignLa lassignLa_br \ +madd mdiv mprod msub mat_binary_op mat_binary_op_br \ +mathprec \ +mcols mcols_br\ +mdiag \ +mdingdong \ +mevsvd mevsvd_br \ +mhilbert \ +mident \ +mlssvd \ +mmult mmult_br \ +mnorms mnorms_br mnormalize mnormalize_br \ +mrange mrange_br \ +mround mround_ij \ +mrows mrows_br\ +msolve msolve_br \ +msum msum_br \ +msvd msvd_br\ +mat_unary_op madjust moffset mscale mat_unary_op_br \ +promote promote_br demote demote_br\ +show show_br \ +transpose transpose_br \ +vdiag vdiag_br \ +vtrim vtrim_br + + global La_namespace + if { $La_namespace != "La" } { + # help advanced users by moving possible compiled code to + # their optionally chosen namespace name. + if { [info commands La::assign_br] == "La::assign_br" } { + rename La:assign_br $La_namespace::assign_br + } + } + +############################# dim #################################### +############################# dim #################################### +############################# dim #################################### +# return the dimension of the argument +# and to some degree verify a proper format +# {} - null +# 0 - scalar +# 1 - vector +# 2 - array +# N - higher +proc dim x { + return [dim_br x] + } + +proc dim_br name_x { + upvar $name_x x + set llen [llength $x] + if { $llen == 0 } {return {}} + if { $llen == 1 && [string is double -strict $x] } { return 0 } + set n [lindex $x 0] + if { $n == 2 } { + set d1 [lindex $x 1] + set d2 [lindex $x 2] + if {[catch { + if { $d1 != int($d1) || $d2 != int($d2) } {error {} } + }] } { error "improper La operand format" } + if { $d2 == 0 } { + if { $llen == [expr {$d1 + 3}] } {return 1 } ;# 2 n 0 + error "improper length of vector" + } + if { $llen == [expr {$d1 * $d2 + 3}] } { return 2 } + error "improper length of matrix" + } + if { $n >= 3 && [string is integer -strict $n] } { + # TBD verify llen as prod of dims + dims + 1 + return $n + } + error "improper La operand format" + } + + + +############################ dotprod ################################# +############################ dotprod ################################# +############################ dotprod ################################# +# +# Dot Product = sum over i, Ai * Bi +# +# Can work columns or rows in matrices +# because indexing increment is passed. +# Default arguments work with vectors, Nx1, or 1XN matrices +# +proc dotprod {a b {N {}} {a0index 3} {b0index 3} {ainc 1} {binc 1}} { + return [dotprod_br a b $N $a0index $b0index $ainc $binc] + } + +# +# always returns a scalar result +# +proc dotprod_br {a_in b_in {N {}} {a0index 3} {b0index 3} {ainc 1} {binc 1}} { + upvar $a_in a + upvar $b_in b + set z 0.0 + if { $N == {}} { ;# default length assumes vector, 1xN, or Nx1 first argument + set len [llength $a] + set N [expr {$len-3}] + # and then see if that matches one of b's dimensions + if { [lindex $b 1] != $N && [lindex $b 2] != $N } { + error "Assumed length does not seem proper." + } + } + for {set i 0} {$i < $N} {incr i} { + set aval [lindex $a $a0index] + set bval [lindex $b $b0index] + set z [expr { $z + $aval * $bval}] + incr a0index $ainc + incr b0index $binc + } + return $z + } + + +########################### join_cols ################################ +########################### join_cols ################################ +########################### join_cols ################################ +# +# combine vector or matrix args as added columns +# +proc join_cols {a b} { + set dimA [dim_br a] + set dimB [dim_br b] + if { $dimA > 2 || $dimB > 2 } { + error "cannot handle > 2D data" + } + if { $dimA < 2 } {promote_br a} + if { $dimB < 2 } {promote_br b} + join_cols_br a b + return $a + } + + +proc join_cols_br {a_in b_in {c_out {}}} { + upvar $a_in a + upvar $b_in b + if { $c_out == {}} {set c_out $a_in } + upvar $c_out c + # have to match row counts + set rows [lindex $a 1] + set acols [lindex $a 2] + set bcols [lindex $b 2] + if { [lindex $b 1] != $rows } { + error "cannot append columns with inequal rows A\[$rows,\] + B\[[lindex $b 1],\]" + } + set total_cols [expr {$acols + $bcols}] + set result [list 2 $rows $total_cols] + for {set row 0} {$row < $rows} {incr row} { + for {set col 0} {$col < $acols} {incr col} { + set i_a [expr {3 + $row * $acols + $col}] + lappend result [lindex $a $i_a] + } + for {set col 0} {$col < $bcols} {incr col} { + set i_b [expr {3 + $row * $bcols + $col}] + lappend result [lindex $b $i_b] + } + } + set c $result + return $c_out + } + +########################### join_rows ################################ +########################### join_rows ################################ +########################### join_rows ################################ +# +# combine vector or matrix args as added rows +# +proc join_rows {a b} { + set dimA [dim_br a] + set dimB [dim_br b] + if { $dimA > 2 || $dimB > 2 } { + error "cannot handle > 2D data" + } + if { $dimA < 2 } {promote_br a} + if { $dimB < 2 } {promote_br b} + join_rows_br a b + return $a + } + + +proc join_rows_br {a_in b_in {c_out {}}} { + upvar $a_in a + upvar $b_in b + if { $c_out == {}} {set c_out $a_in } + upvar $c_out c + # have to match col counts + set arows [lindex $a 1] + set brows [lindex $b 1] + set acols [lindex $a 2] + set bcols [lindex $b 2] + if { $acols != $bcols } { + error "cannot append rows with inequal columns A\[,$acols\] + B\[,$bcols\]" + } + set total_rows [expr {$arows + $brows}] + set c [concat 2 $total_rows $acols [lrange $a 3 end] [lrange $b 3 end]] + return $c_out + } + +########################### lassignLa ################################## +########################### lassignLa ################################## +########################### lassignLa ################################## + +# +# replace a single value in a Tcl list +# This call assumes the caller knows the desired position in the list. +# +proc lassignLa {x index value} { + return [lreplace $x $index $index $value] +} + + +# list element replace by value - this can be done in C code +# without copying the list so it can be efficient, especially since Tcl +# will maintain an internal representation of a list already +# split into elements. +# lreplace forces you to copy the whole list. +if {0} { ;# $tcl_version < 8.4 + # the mysterious K procedure improves efficiency because of how list + # manipulation is performed (thank you Kevin Kenny & Donal Fellows) + proc K {x y} {set x} + + proc lassignLa_br {listname index value} { + upvar $listname thelist + set thelist [lreplace [K $thelist [set thelist {}]] $index $index $value] + return $listname + } +} else { + # lset is available with Tcl 8.4+ + # turns out this is 3X slower! + proc lassignLa_br {listname index value} { + uplevel ::lset $listname $index $value + return $listname + } +} + + +################ madd,mdiv, mprod, msub, mat_binary_op ############## +################ madd,mdiv, mprod, msub, mat_binary_op ############## +################ madd,mdiv, mprod, msub, mat_binary_op ############## +# +# vector/matrix binary operations on elements like addition +# +proc mat_binary_op {a b op} { + set dimA [dim_br a] + set dimB [dim_br b] + if { $dimA > 2 || $dimB > 2 } { + error "cannot handle > 2D data" + } + if { $dimA < 1 || $dimA < $dimB } {promote_br a} + if { $dimB < 1 || $dimB < $dimA } {promote_br b} + mat_binary_op_br a b $op c + return $c + } + +# matrix addition and similar +proc madd {a b} { return [mat_binary_op $a $b +] } +proc mdiv {a b} { return [mat_binary_op $a $b {*1.0/}] } +proc mprod {a b} { return [mat_binary_op $a $b *] } +proc msub {a b} { return [mat_binary_op $a $b -] } + + +# perform a binary operation using corresponding elements of +# two matrices or vectors - op is commonly +, -, *, / +# but you can also add formula constructions as in +# " *2.0 + 3.4*" +# Be careful with division since Tcl tries to perform integer +# division. You may want to use "*1.0 / " to insure floating point math +# +proc mat_binary_op_br {a_in b_in op {c_out {}}} { + upvar $a_in a + upvar $b_in b + if { $c_out == {}} { set c_out $a_in } + upvar $c_out c + set NR [lindex $a 1] + set MC [lindex $a 2] + set NR2 [lindex $b 1] + set MC2 [lindex $b 2] + if { $NR != $NR2 || $MC != $MC2 } { + error "arguments are not conformable A\[$NR,$MC\] vs B\[$NR2,$MC2\]" + } + set result [list 2 $NR $MC] + set imax [llength $a] + for {set i 3} {$i < $imax} {incr i} { + set ai [lindex $a $i] + set bi [lindex $b $i] + lappend result [expr $ai $op $bi] + } + set c $result + return $c_out + } + + +########################## mathprec ################################ +########################## mathprec ################################ +########################## mathprec ################################ +# +# test math precision + +# what is the smallest number such that +# 1+epsilon > 1 +# does the machine truncate or round-off +# what is the radix and number of digits +# +# Intel PC (IEEE 8 byte floating point) results: +# radix=2.0 digits=53 epsilon=2.22044604925e-016 method=truncation +# +proc mathprec {{puts puts}} { + for {set test 1.0} {$test + 1.0 != $test} {set test [expr {$test*2.0}]} { } + for {set diff 1.0} {$test + $diff == $test} {set diff [expr {$diff + 1.0}]} { } + set radix [expr {($test + $diff) - $test}] + if { $diff < $radix } { set method rounding } else {set method truncation } + set digits 0 + for {set test 1.0} {$test + 1 != $test} {set test [expr {$test*$radix}]} { + incr digits + } + set epsilon [expr {pow($radix,(1.0-$digits))}] + if {$puts != {}} { + $puts "radix=$radix digits=$digits epsilon=$epsilon method=$method" + } + return $epsilon + } + +############################ mcols ################################# +############################ mcols ################################# +############################ mcols ################################# +# +# you can use mcols and mrows to determine the number of columns or +# rows and keep your code isolated from the details of the +# data representation +# +proc mcols m { return [lindex $m 2] } +proc mcols_br name_m { upvar $name_m m ; return [lindex $m 2] } +proc mrows m { return [lindex $m 1] } +proc mrows_br name_m { upvar $name_m m ; return [lindex $m 1] } + +############################ mdiag ################################# +############################ mdiag ################################# +############################ mdiag ################################# +# +# create a diagonal matrix from a vector, an Nx1 matrix, or a 1XN matrix +# +proc mdiag v { + set dim [dim_br v] + if { $dim == 2 } { demote_br v } + if { [lindex $v 2] != 0 } { error "expected vector argument" } + set N [lindex $v 1] + set result [list 2 $N $N] + for {set row 0} {$row < $N} {incr row} { + for {set col 0} {$col < $N} {incr col} { + set iv [expr {3 + $row}] + if { $row == $col } { lappend result [lindex $v $iv] }\ + else { lappend result 0 } + } + } + return $result + } + +############################ mdingdong ############################## +############################ mdingdong ############################## +############################ mdingdong ############################## +# +# create the Ding Dong test matrix, a Cauchy matrix that is +# represented inexactly in the machine, but very stable for +# inversion by elimination methods +# +proc mdingdong N { + if { ![string is integer -strict $N] } { error "improper size \"$N\"" } + set result [list 2 $N $N] + for {set row 0} {$row < $N} {incr row} { + for {set col 0} {$col < $N} {incr col} { + lappend result [expr {0.5/($N - $col - $row - 0.5)}] + } + } + return $result + } + +############################ mevsvd ################################# +############################ mevsvd ################################# +############################ mevsvd ################################# +# +# eigenvectors/eigenvalues of a real symmetric matrix by +# singular value decomposition +# +proc mevsvd {A {eps 2.3e-16}} { + mevsvd_br A evals +puts "eigenvalues=[show $evals]" + return $A + } + + +# the eigenvectors are returned in place as the columns of A +# the eigenvalues are returned as a vector + +proc mevsvd_br {A_in_out evals_out {eps 2.3e-16}} { + upvar $A_in_out A + upvar $evals_out evals + set n [lindex $A 1] + if { [lindex $A 2] != $n } { error "expecting square matrix" } + for {set i 0} {$i < $n} {incr i} { + set ii [expr {3 + $i*$n + $i}] + set v [lindex $A $ii] + for {set j 0} {$j < $n} {incr j} { + if { $i != $j } { + set ij [expr {3 + $i*$n + $j}] + set Aij [lindex $A $ij] + set v [expr {$v - abs($Aij)}] + } + } + if { ![info exists h] } { set h $v }\ + elseif { $v < $h } { set h $v } + } + # h is lower bound on Gershgorin region + if { $h <= $eps } { + set h [expr {$h - sqrt($eps)}] + # try to make smallest eigenvalue positive and not too small + for {set i 0} {$i < $n} {incr i} { + set ii [expr {3 + $i*$n + $i}] + set Aii [lindex $A $ii] + lassignLa_br A $ii [expr {$Aii - $h}] + } + }\ + else { + set h 0.0 + } + + set count 0 + # top of the iteration + for {set isweep 0} {$isweep < 30 && $count < $n*($n-1)/2} {incr isweep} { + set count 0 ;# count of rotations in a sweep + + for {set j 0} {$j < [expr {$n-1}]} {incr j} { + for {set k [expr {$j+1}]} {$k < $n} {incr k} { + set p [set q [set r 0.0]] + for {set i 0} {$i < $n} {incr i} { + set ij [expr {3+$i*$n+$j}] + set ik [expr {3+$i*$n+$k}] + set Aij [lindex $A $ij] + set Aik [lindex $A $ik] + set p [expr {$p + $Aij*$Aik}] + set q [expr {$q + $Aij*$Aij}] + set r [expr {$r + $Aik*$Aik}] + } + if { 1.0 >= 1.0 + abs($p/sqrt($q*$r)) } { + if { $q >= $r } { + incr count + # no rotation needed + continue + } + } + set q [expr {$q-$r}] + set v [expr {sqrt(4.0*$p*$p + $q*$q)}] + if { $v == 0.0 } continue + if { $q >= 0.0 } { + set c [expr {sqrt(($v+$q)/(2.0*$v))}] + set s [expr {$p/($v*$c)}] + }\ + else { + set s [expr {sqrt(($v-$q)/(2.0*$v))}] + if { $p < 0.0 } { set s [expr {0.0-$s}] } + set c [expr {$p/($v*$s)}] + } + # rotation + for {set i 0} {$i < $n} {incr i} { + set ij [expr {3+$i*$n+$j}] + set ik [expr {3+$i*$n+$k}] + set Aij [lindex $A $ij] + set Aik [lindex $A $ik] + lassignLa_br A $ij [expr {$Aij*$c + $Aik*$s}] + lassignLa_br A $ik [expr {$Aik*$c - $Aij*$s}] + } + } ;# k + } ;# j + #puts "pass=$isweep skipped rotations=$count" + } ;# isweep + + # now columns are orthogonal, rescale + # and flip signs if all negative or zero + set evals [list 2 $n 0] + for {set j 0} {$j < $n} {incr j} { + set s 0.0 + set notpositive 0 + for {set i 0} {$i < $n} {incr i} { + set ij [expr {3+$i*$n+$j}] + set Aij [lindex $A $ij] + if { $Aij <= 0.0 } { incr notpositive } + set s [expr {$s + $Aij*$Aij}] + } + set s [expr {sqrt($s)}] + if { $notpositive == $n } { set sf [expr {0.0-$s}] }\ + else { set sf $s } + for {set i 0} {$i < $n} {incr i} { + set ij [expr {3+$i*$n+$j}] + set Aij [lindex $A $ij] + lassignLa_br A $ij [expr {$Aij/$sf}] + } + lappend evals [expr {$s+$h}] + } + return $A_in_out + } + +############################ mhilbert ############################### +############################ mhilbert ############################### +############################ mhilbert ############################### +# +# create the Hilbert test matrix which is notorious for being +# ill conditioned for eigenvector/eigenvalue solutions +# +proc mhilbert N { + if { ![string is integer -strict $N] } { error "improper size \"$N\"" } + set result [list 2 $N $N] + for {set row 0} {$row < $N} {incr row} { + for {set col 0} {$col < $N} {incr col} { + lappend result [expr {1.0/($col + $row +1.0)}] + } + } + return $result + } + +############################ mident ################################# +############################ mident ################################# +############################ mident ################################# +# +# create an identity matrix of order N +# +proc mident N { + if { ![string is integer -strict $N] } { error "improper size \"$N\"" } + set result [list 2 $N $N] + for {set row 0} {$row < $N} {incr row} { + for {set col 0} {$col < $N} {incr col} { + if { $row == $col } { lappend result 1 }\ + else { lappend result 0 } + } + } + return $result + } + +############################### mlssvd ######################### +############################### mlssvd ######################### +############################### mlssvd ######################### +# +# solve the over-determined linear equations in the least squares +# sense using SVD +# +# Ax ~ y +# y[m] is dependent variable such as a measured outcome +# x[n] vector of independent variables +# A[m,n] - each row is a set dependent variable values, +# the first column is usually 1 to allow for a constant +# in the regression +# +# q is the minimum singular value, lessor values are treated as 0 +# +proc mlssvd {A y {q 0.0} {puts puts} {epsilon 2.3e-16}} { + # A[m,n] + set m [lindex $A 1] + set n [lindex $A 2] + promote_br y ;# now expect y[m,1] + if { [lindex $y 1] != $m } { + error "cannot conform A\[$m,$n\]*X\[$n] = y\[[lindex $y 1],[lindex $y 2]]" + } + msvd_br A S V + # now A has been transformed to U[m,n] + # S[n], V[n,n] + if { $puts != {}} { + $puts singular\ values=[show $S %.6g] + } + set tol [expr {$epsilon * $epsilon * $n * $n}] + # form Utrans*y into g + set g [list 2 $n 0] + for {set j 0} {$j < $n} {incr j} { + set s 0.0 + for {set i 0} {$i < $m} {incr i} { + set ij [expr {3 + $i*$n + $j}] + set Aij [lindex $A $ij] + set yi [lindex $y [expr {3 + $i}]] + set s [expr {$s + $Aij*$yi}] + } + lappend g $s ;# g[j] = $s + } + # form VS+g = VS+Utrans*g + set x [list 2 $n 0] + for {set j 0} {$j < $n} {incr j} { + set s 0.0 + for {set i 0} {$i < $n} {incr i} { + set iindex [expr {$i+3}] + set zi [lindex $S $iindex] + if { $zi > $q } { + set ji [expr {3 + $j*$n+$i}] + set Vji [lindex $V $ji] + set gi [lindex $g $iindex] + set s [expr {$s + $Vji*$gi/$zi}] + } + } + lappend x $s + } + return $x + } + +############################ mmult ################################## +############################ mmult ################################## +############################ mmult ################################## +# +# matrix multiplication +# +# A[p,q] x B[q,r] = C[p,r] +# Cij = dot product A[i,] row x B[,j] col +# Cij = sum Aik * Bkj ; k=0..q-1 +# +# rules: Vector arguments are always promoted to Nx1 arrays. +# So chances are if you are using one as a left operand +# you probably intend to use the transpose of it (1xN), +# which is easily done using the transpose procedure. +# +proc mmult {a b} { + set dimA [dim_br a] + set dimB [dim_br b] + if { $dimA > 2 || $dimB > 2 } { + error "cannot handle > 2D data" + } + if { $dimA < 2 } {promote_br a} + if { $dimB < 2 } {promote_br b} + mmult_br a b c + return $c + } + +# caller is responsible to provide 2D arguments +proc mmult_br {a_in b_in c_out} { + upvar $a_in a + upvar $b_in b + upvar $c_out c + set p [lindex $a 1] + set q [lindex $a 2] + set q2 [lindex $b 1] + set r [lindex $b 2] + if { $q2 != $q } { + error "matrices are not conformable A\[$p,$q\] x B\[$q2,$r\]" + } + # Cij = dot product A[i,] row x B[,j] col + set c [list 2 $p $r] + set a_base 3 + set ainc 1 + set binc $r + for {set row 0} {$row < $p} {incr row} { + set b_base 3 + for {set col 0} {$col < $r} {incr col} { + lappend c [dotprod_br a b $q $a_base $b_base $ainc $binc] + incr b_base 1 + } + incr a_base $q + } + + return $c_out + } + +######################## mnorms, mnormalize ########################## +######################## mnorms, mnormalize ########################## +######################## mnorms, mnormalize ########################## +# +# compute the means and sigmas of each column +# +# +proc mnorms a { + mnorms_br a means sigmas + return [transpose [join_cols $means $sigmas]] + } + +proc mnorms_br {a_in means_out sigmas_out} { + upvar $a_in a + if { [dim_br a] != 2 } { error "expecting matrix" } + set NR [lindex $a 1] + if { $NR < 2 } { error "insufficient rows to calculate sigmas" } + set NC [lindex $a 2] + # vector results are returned + set means [set sigmas [list 2 $NC 0]] + for {set i 0} {$i < $NC} {incr i} { + mrange_br a colvect $i $i + # mean, stddev are Hume C-code Tcl commands + set coldata [lrange $colvect 3 end] + lappend means [mean $coldata] + lappend sigmas [stddev $coldata] + } + upvar $means_out m + set m $means + upvar $sigmas_out s + set s $sigmas + return [list $means_out $sigmas_out] + } + +# +# normalize each column by subtracting the corresponding mean and then +# dividing by the corresponding sigma +# + +proc mnormalize {a means sigmas} { + if { [dim_br a] == 1 } { promote_br a } + mnormalize_br a means sigmas + return $a + } + +# +# the code will work with means and sigmas as vectors, Nx1, or 1xN matrices +# +proc mnormalize_br {a_in means_in sigmas_in {c_out {}}} { + upvar $a_in a + upvar $means_in means + upvar $sigmas_in sigmas + if { [dim_br a] != 2 } { error "expecting matrix" } + set NR [lindex $a 1] + set NC [lindex $a 2] + set nm [llength $means] + set ns [llength $sigmas] + if { $nm != $ns || $nm != [expr 3+$NC] } { + error "non-conformable arguments a[$NR,$NC] means[expr {$nm-3}] sigmas[expr {$ns-3}]" + } + set result [list 2 $NR $NC] + for {set row 0} {$row < $NR} {incr row} { + for {set i 0} {$i < $NC} {incr i} { + set vindex [expr {3 + $i}] + set mean [lindex $means $vindex] + set sigma [lindex $sigmas $vindex] + set aindex [expr {3 + $row*$NC + $i}] + set val [lindex $a $aindex] + lappend result [expr { (0.0 + $val - $mean)/$sigma }] + } + } + if { $c_out == {}} { set c_out $a_in } + upvar $c_out c + set c $result + return $c_out + } + +# the Hume dmh_wish has these in C code, here they are for other shells +if { [info commands mean] != "mean" } { +proc ::mean numlist { + set len [llength $numlist] + if { $len == 0 } { return {}} + if { $len == 1 } { return [expr {double($numlist)}] } + set s 0.0 + for {set i 0} {$i < $len} {incr i} { + set x [lindex $numlist $i] + set s [expr {$s + $x}] + } + return [expr {$s/$len}] + } +} +if { [info commands stddev] != "stddev" } { +proc ::stddev numlist { + set len [llength $numlist] + if { $len < 2 } { return {} } + set sx [set sxx 0.0] + for {set i 0} {$i < $len} {incr i} { + set x [lindex $numlist $i] + set sx [expr {$sx + $x}] + set sxx [expr {$sxx + $x * $x}] + } + # the real world is full of surprises like stdev {2.41 2.41 2.41} + # causing a domain error, so you get strange code + set diff [expr {$sxx - $sx*$sx/$len}] + if { $diff < 0.0 || ($diff + 0.0 != $diff) } {return 0.0 } + return [expr {sqrt($diff/($len-1.0))}] + } +} + +#################### mrange ########################################## +#################### mrange ########################################## +#################### mrange ########################################## +# +# +# return a subset of selected columns, selected rows +# indexing always starts with 0. "end" can be used as an index. +# You are specifying reverse ordering for the selection when start>last. +# +proc mrange {m colstart collast {rowstart 0} {rowlast end}} { + set dim [dim_br m] + if { $dim < 2 } { promote_br m } + mrange_br m m $colstart $collast $rowstart $rowlast + return $m + } + +# +proc mrange_br {a_in c_out colstart collast {rowstart 0} {rowlast end}} { + upvar $a_in a + # matrix is assumed + set NR [lindex $a 1] + set NC [lindex $a 2] + foreach var {colstart collast} { + if { [set $var] == {end} } { set $var [expr {$NC - 1}] ; continue } + set $var [expr int([set $var])] + if { [set $var] < 0 } { error "column index should not be -ve"}\ + elseif { [set $var] >= $NC } { error "column index [set $var] > [expr {$NC - 1}]" } + } + if { $colstart <= $collast } { set colstep 1 } else {set colstep -1} + foreach var {rowstart rowlast} { + if { [set $var] == {end} } { set $var [expr {$NR - 1}] ; continue } + set $var [expr int([set $var])] + if { [set $var] < 0 } { error "row index should not be -ve"}\ + elseif { [set $var] >= $NR } { error "row index [set $var] > [expr {$NR - 1}]" } + } + if { $rowstart <= $rowlast } { set rowstep 1 } else {set rowstep -1} + set newNR [expr {$rowstep*($rowlast - $rowstart) + 1}] + set newNC [expr {$colstep*($collast - $colstart) + 1}] + set result [list 2 $newNR $newNC] + for {set row $rowstart} {1} {incr row $rowstep} { + for {set col $colstart} {1} {incr col $colstep} { + set index [expr {3 + $row*$NC + $col}] + lappend result [lindex $a $index] + if { $col == $collast } break + } + if { $row == $rowlast } break + } + upvar $c_out c + set c $result + return $c_out + } + +############################ msolve ################################# +############################ msolve ################################# +############################ msolve ################################# + +# +# solve the matrix problem Ax = p for x, where p may be multiple columns +# +# when p is the identity matrix, the solution x, is the inverse of A +# + + +proc msolve {A p {tolerance 2.3e-16}} { + promote_br p ;# a vector becomes Nx1 + + set N [lindex $A 2] + # paste the rhs columns on the end + join_cols_br A p + + msolve_br A $tolerance + + # now peel off the solution which replaced the rhs columns + mrange_br A x $N end + return $x + } + +# +# Gauss elimination with partial pivoting +# +proc msolve_br {A_in {tolerance 2.3e-16}} { + upvar $A_in A + set n [lindex $A 1] + set NC [lindex $A 2] + set p [expr {$NC - $n}] + set Det 1.0 + for {set j 0} {$j <= [expr {$n - 2}]} {incr j} { + set j_plus1 [expr {$j + 1}] + set jj [expr {3 + $j * $NC + $j}] + set Ajj [lindex $A $jj] + set s [expr {abs($Ajj)}] + set k $j + for {set h $j_plus1} {$h < $n} {incr h} { + set Ahj [lindex $A [expr {3+ $h * $NC + $j}]] + set Ahj [expr {abs($Ahj)}] + if { $Ahj > $s } { + set s $Ahj + set k $h + } + } ;# end pivot search + if { $k != $j } { ;# row interchange + for {set i $j} {$i < $NC} {incr i} { + set ki [expr {3+$k*$NC + $i}] + set ji [expr {3+$j*$NC + $i}] + set Aki [lindex $A $ki] + set Aji [lindex $A $ji] + lassignLa_br A $ki $Aji + lassignLa_br A $ji $Aki + } + set Det [expr {0.0-$Det}] + } + set Ajj [lindex $A $jj] + set Det [expr {$Det * $Ajj}] + if { abs($Ajj) <= $tolerance } { + error "matrix is computationally singular at a tolerance of $tolerance" + } + for {set k $j_plus1} {$k < $n} {incr k} { + set kj [expr {3+$k*$NC + $j}] + set Akj [lindex $A $kj] + set Akj [expr {double($Akj)/$Ajj}] + lassignLa_br A $kj $Akj ;# to form multiplier mkj + for {set i $j_plus1} {$i < $NC} {incr i} { + set ki [expr {3+$k*$NC + $i}] + set Aki [lindex $A $ki] + set ji [expr {3+$j*$NC + $i}] + set Aji [lindex $A $ji] + lassignLa_br A $ki [expr {$Aki - double($Akj)*$Aji}] + } + } ;# k + } ;# j + set nn [expr {3+($n-1)*$NC+($n-1)}] ;# n-1,n-1 + set Ann [lindex $A $nn] + set Det [expr {$Det * $Ann}] + if { abs($Ann) < $tolerance } { + error "matrix is computationally singular at a tolerance of $tolerance" + } + # factoring completed + # mij stored in i>j + # Rij stored in i<=j < n + # fij stored in j= n ... n+p-1 + # Back substitution for Rx=f + for {set i $n} {$i < $NC} {incr i} { + set ni [expr {3+($n-1)*$NC + $i}] ;# really n-1,i + set Ani [lindex $A $ni] + set Ani [expr {double($Ani)/$Ann}] + lassignLa_br A $ni $Ani + for {set j [expr {$n-2}]} {$j >= 0} {incr j -1} { + set ji [expr {3+$j*$NC+$i}] + set s [lindex $A $ji] + for {set k [expr {$j+1}]} {$k < $n} {incr k} { + set jk [expr {3+$j*$NC+$k}] + set Ajk [lindex $A $jk] + set ki [expr {3+$k*$NC+$i}] + set Aki [lindex $A $ki] + set s [expr {$s - $Ajk*$Aki}] + } + set jj [expr {3+$j*$NC+$j}] + set Ajj [lindex $A $jj] + lassignLa_br A $ji [expr {double($s)/$Ajj}] + } + } + return $A_in + } + + +########################## msum ##################################### +########################## msum ##################################### +########################## msum ##################################### +# +# compute the sums of each column, return a vector or scalar +# call twice to get sum of columns and rows (set total [msum [msum $a]]) +# +# a is a vector or matrix +# result is a vector or scalar +proc msum a { + if { [dim_br a] == 1 } { promote_br a } + msum_br a sums + # convert vector to scalar if only 1 column + return [demote $sums] + } + +# a_in is matrix +# sums_out writes a vector +proc msum_br {a_in sums_out} { + upvar $a_in a + if { [dim_br a] != 2 } { error "expecting matrix" } + set NR [lindex $a 1] + set NC [lindex $a 2] + # vector results are returned + set sums [list 2 $NC 0] + for {set j 0} {$j < $NC} {incr j} { + set sum 0.0 + set index [expr {3 + $j}] + for {set i 0} {$i < $NR} {incr i} { + set Aij [lindex $a $index] + set sum [expr {$sum + $Aij}] + incr index $NC + } + lappend sums $sum + } + upvar $sums_out m + set m $sums + return $sums_out + } + + +################################ msvd ########################## +################################ msvd ########################## +################################ msvd ########################## +# +# Singular Value Decomposition +# decompose matrix A into (U)(S)(Vtrans) where +# A[m,n] is the original matrix +# U[m,n] has orthogonal columns (Ut)(U) = (1(k) +# and multiplies to an identity matrix ... +# supplemented with zeroes if needed 0(n-k)) +# V[n,n] is orthogonal (V)(Vtran) = [mident $n] +# the eigenvectors representing the principal components +# S is diagonal with the positive singular values of A +# +# when you have V,S,U +# A[m,n]V[n,n] = B[m,n] is a transformation of A to orthogonal columns, B +# B[m,n] = U[m,n]S[n,n] +# square S and divide by (m-1) to get PCA eigenvalues +# +proc msvd A { + msvd_br A S V +puts U:\n[show $A %12.4f] +puts \nS:\n[show $S %12.4f] +puts \nV:\n[show $V %12.4f] + return $V + } +# +proc msvd_br {A_in_U_out S_out V_out {epsilon 2.3e-16}} { + upvar $A_in_U_out A + set m [lindex $A 1] + set n [lindex $A 2] + set tolerance [expr {$epsilon * $epsilon* $m * $n}] + upvar $V_out V + set V [mident $n] + upvar $S_out z + + # top of the iteration + set count 1 + for {set isweep 0} {$isweep < 30 && $count > 0} {incr isweep} { + set count [expr {$n*($n-1)/2}] ;# count of rotations in a sweep + for {set j 0} {$j < [expr {$n-1}]} {incr j} { + for {set k [expr {$j+1}]} {$k < $n} {incr k} { + set p [set q [set r 0.0]] + for {set i 0} {$i < $m} {incr i} { + set ij [expr {3+$i*$n+$j}] + set ik [expr {3+$i*$n+$k}] + set Aij [lindex $A $ij] + set Aik [lindex $A $ik] + set p [expr {$p + $Aij*$Aik}] + set q [expr {$q + $Aij*$Aij}] + set r [expr {$r + $Aik*$Aik}] + } + if { $q < $r } { + set c 0.0 + set s 1.0 + }\ + elseif { $q * $r == 0.0 } { ;# underflow of small elements + incr count -1 + continue + }\ + elseif { ($p*$p)/($q*$r) < $tolerance } { ;# cols j,k are orthogonal + incr count -1 + continue + }\ + else { + set q [expr {$q-$r}] + set v [expr {sqrt(4.0*$p*$p + $q*$q)}] + set c [expr {sqrt(($v+$q)/(2.0*$v))}] + set s [expr {$p/($v*$c)}] + # s == sine of rotation angle, c == cosine + } + # rotation of A + for {set i 0} {$i < $m} {incr i} { + set ij [expr {3+$i*$n+$j}] + set ik [expr {3+$i*$n+$k}] + set Aij [lindex $A $ij] + set Aik [lindex $A $ik] + lassignLa_br A $ij [expr {$Aij*$c + $Aik*$s}] + lassignLa_br A $ik [expr {$Aik*$c - $Aij*$s}] + } + # rotation of V + for {set i 0} {$i < $n} {incr i} { + set ij [expr {3+$i*$n+$j}] + set ik [expr {3+$i*$n+$k}] + set Vij [lindex $V $ij] + set Vik [lindex $V $ik] + lassignLa_br V $ij [expr {$Vij*$c + $Vik*$s}] + lassignLa_br V $ik [expr {$Vik*$c - $Vij*$s}] + } + } ;# k + } ;# j + #puts "pass=$isweep skipped rotations=$count" + } ;# isweep + + set z [list 2 $n 0] + for {set j 0} {$j < $n} {incr j} { + set q 0.0 + for {set i 0} {$i < $m} {incr i} { + set ij [expr {3+$i*$n+$j}] + set Aij [lindex $A $ij] + set q [expr {$q+$Aij*$Aij}] + } + set q [expr {sqrt($q)}] + lappend z $q + if { $q >= $tolerance } { + for {set i 0} {$i < $m} {incr i} { + set ij [expr {3+$i*$n+$j}] + set Aij [lindex $A $ij] + lassignLa_br A $ij [expr {$Aij/$q}] + } + } + } ;# j + return [list $A_in_U_out $S_out $V_out] + } + + +######### mat_unary_op, madjust, moffset, mscale #################### +######### mat_unary_op, madjust, moffset, mscale #################### +######### mat_unary_op, madjust, moffset, mscale #################### +# +# matrix unary operations like scaling +# +proc mat_unary_op {a op} { + set dimA [dim_br a] + if { $dimA > 2} { + error "cannot handle > 2D data" + } + if { $dimA < 1 } { promote_br a } + mat_unary_op_br a $op c + return $c + } + +proc moffset {a delta} {return [mat_unary_op $a "expr $delta +"] } +proc mscale {a delta} { return [mat_unary_op $a "expr $delta *"] } +proc madjust {a scale offset} { return [mat_unary_op $a "expr $offset + $scale *"] } +proc mround_ij {eps aij} { ;# round off a number if close to an integer + set delta [expr {abs($aij - round($aij))}] + if { $delta <= $eps || $delta <= abs($eps * $aij) } { return [expr {double(round($aij))}] } + return $aij + } +proc mround {a {eps 1.0e-8}} { return [mat_unary_op $a [list mround_ij $eps]] } + +# perform a unary operation on the elements of a vector or matrix +# the value is lappended to your op and the result is executed +# "expr 0.5 *" divides by 2 +# you can use a procedure for more complex formulae +# +proc mat_unary_op_br {a_in op {c_out {}}} { + upvar $a_in a + if { $c_out == {}} { set c_out $a_in } + upvar $c_out c + set NR [lindex $a 1] + set MC [lindex $a 2] + set result [list 2 $NR $MC] + set imax [llength $a] + for {set i 3} {$i < $imax} {incr i} { + lappend result [eval [concat $op [lindex $a $i]]] + } + set c $result ;# only changes a_in,c if no error + return $c_out + } +########################### promote ################################# +########################### promote ################################# +########################### promote ################################# +# +# promote a scalar or vector to an array +# vec[N] promoted to Nx1 array +# +proc promote x { + promote_br x + return $x + } + +proc promote_br {name_in {name_out {}}} { + upvar $name_in x + if { $name_out == {}} {set name_out $name_in } + upvar $name_out y + set dim [dim_br x] + if { $dim == 0 } { + set y [list 2 1 1 $x] + return $name_out + } + if { $dim == 1 } { ;# 2 n 0 v1 v2 ... + set n [lindex $x 1] + if { $name_out != $name_in } { set y [lreplace $x 2 2 1] }\ + else { lassignLa_br y 2 1 } + return $name_out + } + if { $name_out != $name_in } {set y $x} + return $name_out + } + + +# demote an Nx1 or 1XN matrix to a vect[n] +# demote a vect[1] to a scalar +# call twice to demote 1x1 matrix to scalar +proc demote x { + demote_br x + return $x + } + +proc demote_br {a_in {c_out {}}} { + upvar $a_in a + if { $c_out == {}} { set c_out $a_in } + upvar $c_out c + if { $c_out != $a_in } { set c $a } + # demote 1XN or Nx1 to vector + set dim [dim_br a] + if { $dim == 0 || $dim > 2 } {return $c_out } + set nr [lindex $a 1] + if { $dim == 2 } { + set mc [lindex $a 2] + if { $mc == 1 } { + set c [lreplace $a 2 2 0] + }\ + elseif { $nr == 1 } { + set c [lreplace $a 1 2 $mc 0] + } + return $c_out + } + # demote vector1 to scalar + if { $dim == 1 } { + if { $nr == 1 } { + set c [lindex $a 3] + } + } + return $c_out + } + +############################ show ################################### +############################ show ################################### +############################ show ################################### + +# +# Return a formatted string representation for an lalf operand +# the format argument is optionally used per-element +# with the Tcl format command. Examples: +# %.4f gives 4 decimal digit fixed point +# %12.4e provides 4 decimal digit exponential format with a fixed +# width +# +# The col_join argument can be used to separate values with commas +# and tabs, etc. Similarly the row_join allows you to define the +# row separation string. +# +proc show {x {format {}} {col_join { }} {row_join \n}} { + show_br x x $format $col_join $row_join + return $x + } + +proc show_br {name_in {name_out {}} {format {}} {col_join { }} {row_join \n}} { + upvar $name_in x + if { $name_out == {}} { set name_out $name_in } + upvar $name_out result + set dim [dim_br x] + if { $dim == 0 } { + if { $format != {} } {set result [format $format $x] }\ + else { set result $x } + return $name_out} + if { $dim == 1 } { ;# 2 n 0 v1 v2 ... + if { $format == {}} { + set result [join [lrange $x 3 end] $col_join] + }\ + else { + set temp [list] + foreach item [lrange $x 3 end] { + lappend temp [format $format $item] + } + set result [join $temp $col_join] + } + return $name_out + } + if { $dim == 2} { + if { $format != {}} { + mat_unary_op_br x [list format $format] result + } \ + else { + set result $x + } + set rows [list] + set NR [lindex $result 1] + set MC [lindex $result 2] + for {set row 0} {$row < $NR} {incr row} { + set start [expr {3 + $row*$MC}] + set end [expr {$start + $MC - 1}] + lappend rows [join [lrange $result $start $end] $col_join] + } + set result [join $rows $row_join] + return $name_out + } + + error "> 2-D TBD" + } + + +########################### transpose ############################### +########################### transpose ############################### +########################### transpose ############################### +# +# exchange [i,j] with [j,i] +# +# a vector is promoted to a 1xN array by transpose (1 row x N col) +# +#% set a {2 2 3 01 02 03 11 12 13} +#% show $a +# 01 02 03 +# 11 12 13 +#% transpose $a +# 2 3 2 01 11 02 12 03 13 +#% show [transpose $a] +# 01 11 +# 02 12 +# 03 13 +# +proc transpose x { ;# call by value + transpose_br x + return $x + } + +# transpose +# call by reference - default output is to overwrite input data +# +proc transpose_br {name_in {name_out {}}} { + upvar $name_in x + if { $name_out == {}} {set name_out $name_in } + upvar $name_out y + set dim [dim_br x] + if { $dim == 0 } { return } + if { $dim == 1 } { ;# {2 N 0 ...} vector to 1xN + set N [lindex $x 1] + set y [lreplace $x 1 2 1 $N] + return $name_out + } + if { $dim == 2 } { + + set NR [lindex $x 1] + set MC [lindex $x 2] + set result [list 2 $MC $NR] + for {set j 0} {$j < $MC} {incr j} { ;# j == col in source + for { set i 0 } { $i < $NR } { incr i } { ;# i == row in source + lappend result [lindex $x [expr {3 + $i * $MC + $j}]] + } + } + set y $result + return $name_out + } + error "transpose cannot handle >2D data" + } + +########################### vdiag ############################# +########################### vdiag ############################# +########################### vdiag ############################# +# +# Create a vector from the diagonal elements of a matrix +# +proc vdiag m { + promote_br m + vdiag_br m v + return $v + } + +proc vdiag_br {name_m_in name_v_out} { + upvar $name_m_in m + upvar $name_v_out v + set N [lindex $m 1] + set C [lindex $m 2] + if { $C < $N } { set N $C } + set v [list 2 $N 0] + for {set i 0} {$i < $N} {incr i} { + set ij [expr {$i*($C + 1) + 3}] + lappend v [lindex $m $ij] + } + + return $name_v_out + } + + +########################### vtrim ############################# +########################### vtrim ############################# +########################### vtrim ############################# + +# +# for a vector, just return the actual vector elements +# trim away the dimension and size data in the front +# +proc vtrim x { + vtrim_br x + return $x + } + +proc vtrim_br {name_in {name_out {}}} { + upvar $name_in x + if { $name_out == {}} {set name_out $name_in } + upvar $name_out y + set y [lrange $x 3 end] + return $name_out + } + + + + + + + + + + + +} ;# end of namespace diff --git a/lib/la1.0/la.test b/lib/la1.0/la.test new file mode 100755 index 0000000000000000000000000000000000000000..95c4190493d55792310e52140c3b9ec507043ffb --- /dev/null +++ b/lib/la1.0/la.test @@ -0,0 +1,500 @@ +# Part of the Hume La package +# (C)Copyright 2001, Hume Integration Software + +source la.tcl +catch { namespace import La::*} + +if { [info commands dim] == {} } { + puts "Cannot find La package procedures like dim" + return + } + +if {[string compare test [info procs test]] == 1} then \ + {source defs} + +test mathprec {compute precision of math calculations} { + set eps [mathprec] + expr $eps + 1 > 1 +} {1} + +test dim-1 {various good args} { + set list {} + foreach arg { {} {2.3} {2 1 0 2.3} {2 2 2 1.1 1.2 1.3 1.4} {3 1 1 1 1}} { + lappend list [dim $arg] + } + set list +} {{} 0 1 2 3} + +test dim-2 {bad vector arg} { + catch {dim {2 0 3 1 2 3}} txt + set txt +} {improper length of matrix} + +test dim-3 {bad scalar arg} { + catch {dim {2 0 0 3.14}} txt + set txt +} {improper length of vector} + +test dim-4 {bad vector arg} { + catch {dim {2 3 0 1 2 3 4}} txt + set txt +} {improper length of vector} + +test dim-5 {overlong matrix arg} { + catch {dim {2 3 1 1 2 3 4}} txt + set txt +} {improper length of matrix} + +test dim-6 {short matrix arg} { + catch {dim {2 3 1 1 2}} txt + set txt +} {improper length of matrix} + +test dim-7 {incomplete vector/matrix} { + catch {dim {2 1}} txt + set txt +} {improper La operand format} + +test dim-8 {wrong type arg} { + catch {dim {hello world}} txt + set txt +} {improper La operand format} + +test dim-9 {wrong type arg} { + catch {dim {hello}} txt + set txt +} {improper La operand format} + +test dotprod-1 {vxv good} { + dotprod {2 3 0 1 2 3} {2 3 0 1 2 3} +} {14.0} + +test dotprod-2 {v x v bad} { + catch {dotprod {2 3 0 1 2 3} {2 4 0 1 2 3 4}} txt + set txt +} {Assumed length does not seem proper.} + +test dotprod-3 {v x Nx1 good} { + dotprod {2 3 0 1 2 3} {2 3 1 1 2 3} +} {14.0} + +test dotprod-4 {v x 1XN good} { + dotprod {2 3 0 1 2 3} {2 1 3 1 2 3} +} {14.0} + +test dotprod-5 {1xN x Nx1 good} { + dotprod {2 1 3 1 2 3} {2 3 1 1 2 3} +} {14.0} + +test dotprod-6 {1xN x v good} { + dotprod {2 1 3 1 2 3} {2 3 0 1 2 3} +} {14.0} + +test dotprod-7 {1xN x 1xN good} { + dotprod {2 1 3 1 2 3} {2 1 3 1 2 3} +} {14.0} +test dotprod-8 {Nx1 x v good} { + dotprod {2 3 1 1 2 3} {2 3 0 1 2 3} +} {14.0} +test dotprod-9 {Nx1 x 1xN good} { + dotprod {2 3 1 1 2 3} {2 1 3 1 2 3} +} {14.0} +test dotprod-10 {Nx1 x Nx1 good} { + dotprod {2 3 1 1 2 3} {2 3 1 1 2 3} +} {14.0} + +test dotprod-11 {Nx1 x NxN good} { + dotprod {2 3 1 1 2 3} {2 3 3 1 2 3 4 5 6 7 8 9} +} {14.0} + +test dotprod-12 {Nx1 x NxN good, 2nd row} { + dotprod {2 3 1 1 2 3} {2 3 3 1 2 3 0 1 0 7 8 9} 3 3 6 +} {2.0} + +test dotprod-13 {Nx1 x NxN good, 2nd col} { + dotprod {2 3 1 1 2 3} {2 3 3 1 2 3 0 1 0 7 8 9} 3 3 4 1 3 +} {28.0} + +test join_cols-1 {v v good} { + join_cols {2 3 0 1 2 3} {2 3 0 1.1 2.2 3.3} +} {2 3 2 1 1.1 2 2.2 3 3.3} + +test join_cols-2 {v v bad} { + catch {join_cols {2 3 0 1 2 3} {2 4 0 1.1 2.2 3.3 4.4}} txt + set txt +} {cannot append columns with inequal rows A[3,] + B[4,]} + +test join_cols-3 {m v good} { + join_cols {2 3 3 1 2 3 4 5 6 7 8 9} {2 3 0 1.1 2.2 3.3} +} {2 3 4 1 2 3 1.1 4 5 6 2.2 7 8 9 3.3} + +test join_cols-4 {v m good} { + join_cols {2 3 0 1.1 2.2 3.3} {2 3 3 1 2 3 4 5 6 7 8 9} +} {2 3 4 1.1 1 2 3 2.2 4 5 6 3.3 7 8 9} + +test join_cols-3 {m m good} { + join_cols {2 3 3 1 2 3 4 5 6 7 8 9} {2 3 3 1 2 3 4 5 6 7 8 9} +} {2 3 6 1 2 3 1 2 3 4 5 6 4 5 6 7 8 9 7 8 9} + +test join_rows-1 {v v equal} { + join_rows {2 3 0 1 2 3} {2 3 0 1.1 2.2 3.3} +} {2 6 1 1 2 3 1.1 2.2 3.3} + +test join_rows-2 {v v not equal} { + catch {join_rows {2 3 0 1 2 3} {2 2 0 1.1 2.2}} txt + set txt +} {2 5 1 1 2 3 1.1 2.2} + +test join_rows-3 {m v bad} { + catch {join_rows {2 3 3 1 2 3 4 5 6 7 8 9} {2 3 0 1.1 2.2 3.3}} txt + set txt +} {cannot append rows with inequal columns A[,3] + B[,1]} + +test join_rows-4 {m v good} { + catch {join_rows {2 3 3 1 2 3 4 5 6 7 8 9} [transpose {2 3 0 1.1 2.2 3.3}]} txt + set txt +} {2 4 3 1 2 3 4 5 6 7 8 9 1.1 2.2 3.3} + +test join_rows-5 {v m good} { + join_rows [transpose {2 3 0 1.1 2.2 3.3}] {2 3 3 1 2 3 4 5 6 7 8 9} +} {2 4 3 1.1 2.2 3.3 1 2 3 4 5 6 7 8 9} + +test join_cols-6 {m m good} { + join_rows {2 3 3 1 2 3 4 5 6 7 8 9} {2 3 3 1 2 3 4 5 6 7 8 9} +} {2 6 3 1 2 3 4 5 6 7 8 9 1 2 3 4 5 6 7 8 9} + +test lassign_br-1 {good} { + set a {1 2 3} + lassign_br a 1 2.2 + set a +} {1 2.2 3} + +test lassign_br-2 {out-of-range} { + set a {1 2 3} + catch {lassign_br a 3 2.2} txt +} {1} + +test madd-1 {matrix add} { + set a {2 2 2 1 2 3 4} + madd $a $a +} {2 2 2 2 4 6 8} + +test madd-2 {matrix add} { + set a {2 2 2 1 2 3 4} + catch {madd $a {2 1 1 2}} txt + set txt +} {arguments are not conformable A[2,2] vs B[1,1]} + +test mcols-1 {col count} { + mcols {2 1 3 1 2 3} +} {3} + + +test mdiag-1 {v} { + mdiag {2 3 0 1 2 3} +} {2 3 3 1 0 0 0 2 0 0 0 3} + +test mdiag-2 {1xN} { + mdiag {2 1 3 1 2 3} +} {2 3 3 1 0 0 0 2 0 0 0 3} + +test mdiag-3 {Nx1} { + mdiag {2 3 1 1 2 3} +} {2 3 3 1 0 0 0 2 0 0 0 3} + + +test mevsvd-1 {evals of hilbert matrix} { + global mevect + set mevect [mhilbert 4] + mevsvd_br mevect evals + set txt evals=\n[show $evals %.6f]\n + append txt vects=\n[show $mevect %.6f] + set txt +} {evals= +1.500214 0.169141 0.006738 0.000097 +vects= +0.792608 -0.582076 0.179186 -0.029193 +0.451923 0.370502 -0.741918 0.328712 +0.322416 0.509579 0.100228 -0.791411 +0.252161 0.514048 0.638283 0.514553} + +test mevsvd-2 {evect test} { + global mevect + set m [mmult [transpose $mevect] $mevect] + show [mround $m] +} {1.0 0.0 0.0 0.0 +0.0 1.0 0.0 0.0 +0.0 0.0 1.0 0.0 +0.0 0.0 0.0 1.0} + + +test mident-1 {1} { + mident 1 +} {2 1 1 1} + +test mident-2 {3} { + mident 3 +} {2 3 3 1 0 0 0 1 0 0 0 1} + +test mident-3 {3.5} { + catch {mident 3.5} txt + set txt +} {improper size "3.5"} + +test mlssvd {Nash data} { + set A {2 13 5 \ + 1 563 262 461 221\ + 1 658 291 473 222\ + 1 676 294 513 221\ + 1 749 302 516 218\ + 1 834 320 540 217\ + 1 973 350 596 218\ + 1 1079 386 650 218\ + 1 1151 401 676 225\ + 1 1324 446 769 228\ + 1 1499 492 870 230\ + 1 1690 510 907 237\ + 1 1735 534 932 235\ + 1 1778 559 956 236} + set y {2 13 0\ + 305 342 331 339 354 369 378 368 405 438 438 451 485} + set x [mlssvd $A $y 0.0 {}] + show $x %.6f +} {207.782626 -0.046192 1.019387 -0.159823 -0.290376} + +test mmult-1 {good} { + set id3 [mident 3] + mmult $id3 $id3 +} {2 3 3 1.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 1.0} + +test mmult-2 {good} { + set id3 [mident 3] + set m {2 3 3 1 2 3 4 5 6 7 8 9} + mmult $id3 $m +} {2 3 3 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0} + +test mmult-3 {good} { + set id3 [mident 3] + set m {2 3 3 1 2 3 4 5 6 7 8 9} + mmult $m $id3 +} {2 3 3 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0} + +test mmult-4 {m v good} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + mmult $m {2 3 0 1 2 3} +} {2 3 1 14.0 32.0 50.0} + +test mmult-5 {m v mismatch} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + catch {mmult $m {2 4 0 1 2 3 4}} txt + set txt +} {matrices are not conformable A[3,3] x B[4,1]} + +test mmult-6 {v m mismatch} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + catch {mmult {2 3 0 1 2 3} $m} txt + set txt +} {matrices are not conformable A[3,1] x B[3,3]} + +test mmult-7 {vt m match} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + catch {mmult [transpose {2 3 0 1 2 3}] $m} txt + set txt +} {2 1 3 30.0 36.0 42.0} + +test mnorms-1 {id3} { + show [mnorms [mident 3]] %6.4f +} {0.3333 0.3333 0.3333 +0.5774 0.5774 0.5774} + +test mnorms-2 {3x3 mat} { + show [mnorms {2 3 3 1 2 3 4 5 6 7 8 9}] %6.4f +} {4.0000 5.0000 6.0000 +3.0000 3.0000 3.0000} + +test mnorms-3 {2x2 mat} { + show [mnorms {2 2 2 1 2 1 2}] %6.4f +} {1.0000 2.0000 +0.0000 0.0000} + +test mnorms-4 {4x2 mat} { + show [mnorms {2 4 2 1 2 -1 -2 3 5 -5 -3}] %6.4f +} {-0.5000 0.5000 +3.4157 3.6968} + +test mnormalize-1 {4x2 mat} { + set m {2 4 2 1 2 -1 -2 3 5 -5 -3} + mnorms_br m means sigmas + show [mnormalize $m $means $sigmas] %6.4f +} {0.4392 0.4058 +-0.1464 -0.6763 +1.0247 1.2173 +-1.3175 -0.9468} + +test mrange-1 {3x3 mat} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + mrange $m 0 1 +} {2 3 2 1 2 4 5 7 8} + +test mrange-2 {3x3 mat} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + mrange $m 2 2 +} {2 3 1 3 6 9} + +test mrange-3 {3x3 mat index err} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + catch {mrange $m 3 3} txt + set txt +} {column index 3 > 2} + +test mrows-1 {row count} { + mrows {2 1 3 1 2 3} +} {1} + +test msolve-1 {inverse of hilbert matrix} { + global mevect + set mevect [mhilbert 4] + set mevect [msolve $mevect [mident 4]] + show [mround $mevect] +} {16.0 -120.0 240.0 -140.0 +-120.0 1200.0 -2700.0 1680.0 +240.0 -2700.0 6480.0 -4200.0 +-140.0 1680.0 -4200.0 2800.0} + +test msolve-2 {minv test} { + global mevect + set m [mmult [transpose [mhilbert 4]] $mevect] + show [mround $m] +} {1.0 0.0 0.0 0.0 +0.0 1.0 0.0 0.0 +0.0 0.0 1.0 0.0 +0.0 0.0 0.0 1.0} + +test moffset-1 {add const} { + set m {2 3 2 1 2 3 4 5 6} + moffset $m 0 +} {2 3 2 1 2 3 4 5 6} + +test moffset-2 {add const} { + set m {2 3 2 1 2 3 4 5 6} + moffset $m 0.1 +} {2 3 2 1.1 2.1 3.1 4.1 5.1 6.1} + +test mscale-1 {mult const} { + set m {2 3 2 1 2 3 4 5 6} + mscale $m 0.1 +} {2 3 2 0.1 0.2 0.3 0.4 0.5 0.6} + +test madjust-1 {mult const} { + set m {2 3 2 1 2 3 4 5 6} + madjust $m 2.0 0.1 +} {2 3 2 2.1 4.1 6.1 8.1 10.1 12.1} + +test msum-1 {reduce vector} { + msum {2 3 0 1 2 3} +} {6.0} + +test msum-2 {reduce matrix cols} { + msum {2 3 3 1 2 3 4 5 6 7 8 9} +} {2 3 0 12.0 15.0 18.0} + +test msum-3 {reduce matrix to scalar} { + msum [msum {2 3 3 1 2 3 4 5 6 7 8 9}] +} {45.0} + +test promote-1 {scalar} { + promote 3.1415 +} {2 1 1 3.1415} + +test promote-2 {vector} { + promote {2 2 0 1 2} +} {2 2 1 1 2} + +test promote-3 {matrix} { + promote {2 2 1 1 2} +} {2 2 1 1 2} + +test promote-4 {3-D} { + promote {3 1 1 1 3.14} +} {3 1 1 1 3.14} + +test demote-1 scalar { + demote 3.14 +} {3.14} + +test demote-2 {vector 1} { + demote {2 1 0 3.14} +} {3.14} + +test demote-3 {vector 2} { + demote {2 2 0 1 2} +} {2 2 0 1 2} + +test demote-4 {matrix Nx1} { + demote {2 3 1 1 2 3} +} {2 3 0 1 2 3} + +test demote-5 {matrix 1xN} { + demote {2 1 3 1 2 3} +} {2 3 0 1 2 3} + +test demote-6 {garbarge} { + catch {demote {garbage in garbage out}} txt + set txt +} {improper La operand format} + +test show-1 {optional separators} { + set m {2 2 3 1 2 3 4 5 6} + show $m {} , \; +} {1,2,3;4,5,6} + +test show-2 {optional separators, non default format} { + set m {2 2 3 1 2 3 4 5 6} + show $m %g , \; +} {1,2,3;4,5,6} + +test show-3 {optional separators, non-default format} { + set m {2 2 3 1 2 3 4 5 6} + show $m %.1f , \; +} {1.0,2.0,3.0;4.0,5.0,6.0} + +test show-4 {optional separators} { + set v {2 3 0 1 2 3} + show $v {} , +} {1,2,3} + +test show-5 {optional separators, non default format} { + set v {2 3 0 1 2 3} + show $v %g , +} {1,2,3} + +test show-6 {optional separators, non default format} { + set v {2 3 0 1 2 3} + show $v %.1f , +} {1.0,2.0,3.0} + +test transpose-1 vector { + transpose {2 3 0 1 2 3} +} {2 1 3 1 2 3} + +test transpose-2 {nxm mat} { + transpose {2 3 2 1 2 3 4 5 6} +} {2 2 3 1 3 5 2 4 6} + +test vdiag-1 {diag of matrix} { + set m {2 3 3 1 2 3 4 5 6 7 8 9} + vdiag $m +} {2 3 0 1 5 9} + +test vdiag-2 {diag of matrix} { + set m {2 2 3 1 2 3 4 5 6} + vdiag $m +} {2 2 0 1 5} + +test vtrim-1 vect { + vtrim {2 3 0 1 2 3} +} {1 2 3} + + + + diff --git a/lib/la1.0/la_cmds.html b/lib/la1.0/la_cmds.html new file mode 100755 index 0000000000000000000000000000000000000000..4006b5b7d6b8653e968fa724a9732deffd8a7ded --- /dev/null +++ b/lib/la1.0/la_cmds.html @@ -0,0 +1,113 @@ + + +Linear Algebra (La) Package Contents + +

+

Contents

+

+ +

+La Feature Summary +

+

User Guide

+ +Obtaining +

+Installation +

+Runtime Usage +

+Operand Formats +

+Argument Passing +

+Performance +

+PCA Example +

+

Reference

+

+demote +
+dim +
+dotprod +
+join_cols +
+join_rows +
+lassign +
+madd +
+madjust +
+mat_binary_op +
+mat_unary_op +
+mathprec +
+mcols +
+mdiag +
+mdingdong +
+mdiv +
+mevsvd +
+mhilbert +
+mident +
+mlssvd +
+mmult +
+mnorms +
+mnormalize +
+moffset +
+mprod +
+mrange +
+mround +
+mrows +
+mscale +
+msolve +
+msub +
+msum +
+msvd +
+promote +
+show +
+transpose +
+vdiag +
+vtrim + +

+

Author

+
+

References

+
+

License Terms

+
+

Document Version

+ + diff --git a/lib/la1.0/la_download.html b/lib/la1.0/la_download.html new file mode 100755 index 0000000000000000000000000000000000000000..526bdfbd7bdcf3f0ac926f47b4267bf2b5a050dd --- /dev/null +++ b/lib/la1.0/la_download.html @@ -0,0 +1,28 @@ + + +Tcl Linear Algebra (La) Package Download + +

+

Linear Algebra Package Download

+

+La package, version 1.01, for Windows (text files have CR/LF line delimiters) +
+

+

+La package for POSIX (text files have LF line delimiters)
+

+

+

Version History

+Version 1.01 differs from version 1.0 only in the renaming of +Hume Integration Services to Hume Integration Software. +

+

+


+The Hume Integration Software Home Page. + + + diff --git a/lib/la1.0/license.txt b/lib/la1.0/license.txt new file mode 100755 index 0000000000000000000000000000000000000000..5b86a3505815bc6a3bc32acd1b2f3c38bfd03c66 --- /dev/null +++ b/lib/la1.0/license.txt @@ -0,0 +1,44 @@ +The La package software is being distributed under terms and +conditions similar to Tcl/Tk. The author is providing the +package software to the Tcl community as a returned favor for +the value of packages received over the years. + + +The La package software is copyrighted by Hume Integration Software. +The following terms apply to all files associated with the software +unless explicitly disclaimed in individual files. + +The authors hereby grant permission to use, copy, modify, distribute, +and license this software and its documentation for any purpose, provided +that existing copyright notices are retained in all copies and that this +notice is included verbatim in any distributions. No written agreement, +license, or royalty fee is required for any of the authorized uses. +Modifications to this software may be copyrighted by their authors +and need not follow the licensing terms described here, provided that +the new terms are clearly indicated on the first page of each file where +they apply. + +IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY +FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES +ARISING OUT OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY +DERIVATIVES THEREOF, EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE. + +THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, +INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE, AND NON-INFRINGEMENT. THIS SOFTWARE +IS PROVIDED ON AN "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE +NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR +MODIFICATIONS. + +GOVERNMENT USE: If you are acquiring this software on behalf of the +U.S. government, the Government shall have only "Restricted Rights" +in the software and related documentation as defined in the Federal +Acquisition Regulations (FARs) in Clause 52.227.19 (c) (2). If you +are acquiring the software on behalf of the Department of Defense, the +software shall be classified as "Commercial Computer Software" and the +Government shall have only "Restricted Rights" as defined in Clause +252.227-7013 (c) (1) of DFARs. Notwithstanding the foregoing, the +authors grant the U.S. Government and others acting in its behalf +permission to use and distribute the software in accordance with the +terms specified in this license. diff --git a/lib/la1.0/pkgIndex.tcl b/lib/la1.0/pkgIndex.tcl new file mode 100755 index 0000000000000000000000000000000000000000..8314d214b3b74655ea7be80d0e45e0776730e27c --- /dev/null +++ b/lib/la1.0/pkgIndex.tcl @@ -0,0 +1,11 @@ +# Tcl package index file, version 1.1 +# This file is generated by the "pkg_mkIndex" command +# and sourced either when an application starts up or +# by a "package unknown" script. It invokes the +# "package ifneeded" command to set up package-related +# information so that packages will be loaded automatically +# in response to "package require" commands. When this +# script is sourced, the variable $dir must contain the +# full path name of this file's directory. + +package ifneeded La 1.0.1 [list source [file join $dir la.tcl]] diff --git a/lib/la1.0/products.html b/lib/la1.0/products.html new file mode 100755 index 0000000000000000000000000000000000000000..944b97721207ac7ef4d2998693aa2660aed4576f --- /dev/null +++ b/lib/la1.0/products.html @@ -0,0 +1,15 @@ + + +Hume Integration Tcl/Tk/DMH Documentation + + +

+

Hume Integration Software

+

Document Sections

+ +Linear Algebra Package
+

+ +LA Package Download
+ + \ No newline at end of file diff --git a/lib/la1.0/tclIndex b/lib/la1.0/tclIndex new file mode 100755 index 0000000000000000000000000000000000000000..733aabb1db71c325c243131801d2f79835abd06a --- /dev/null +++ b/lib/la1.0/tclIndex @@ -0,0 +1,69 @@ +# Tcl autoload index file, version 2.0 +# This file is generated by the "auto_mkindex" command +# and sourced to set up indexing information for one or +# more commands. Typically each line is a command that +# sets an element in the auto_index array, where the +# element name is the name of a command and the value is +# a script that loads the command. + +set {auto_index(::$La_namespace::dim)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::dim_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::dotprod)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::dotprod_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::join_cols)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::join_cols_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::join_rows)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::join_rows_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::lassign)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mat_binary_op)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::madd)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mdiv)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mprod)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msub)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mat_binary_op_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mathprec)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mcols)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mcols_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mrows)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mrows_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mdiag)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mdingdong)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mevsvd)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mevsvd_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mhilbert)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mident)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mlssvd)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mmult)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mmult_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mnorms)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mnorms_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mnormalize)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mnormalize_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mrange)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mrange_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msolve)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msolve_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msum)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msum_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msvd)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::msvd_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mat_unary_op)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::moffset)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mscale)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::madjust)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mround_ij)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mround)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::mat_unary_op_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::promote)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::promote_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::demote)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::demote_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::show)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::show_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::transpose)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::transpose_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::vdiag)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::vdiag_br)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::vtrim)} [list source [file join $dir la.tcl]] +set {auto_index(::$La_namespace::vtrim_br)} [list source [file join $dir la.tcl]] +set auto_index(nist_pca) [list source [file join $dir NIST.tcl]]