| 1 | # $Header: /usr/cvsroot/la/la.tcl,v 1.2 2001/07/12 23:59:51 hume Exp $ |
|---|
| 2 | # |
|---|
| 3 | # La == Linear Algebra and similar math |
|---|
| 4 | # |
|---|
| 5 | # Part of the La Package |
|---|
| 6 | # (C)Copyright 2001, Hume Integration Services |
|---|
| 7 | # |
|---|
| 8 | # $Log: la.tcl,v $ |
|---|
| 9 | # Revision 1.2 2001/07/12 23:59:51 hume |
|---|
| 10 | # Added K proc and boosted performance bigtime! |
|---|
| 11 | # |
|---|
| 12 | # Revision 1.1.1.1 2001/07/12 15:50:16 hume |
|---|
| 13 | # Initial check-in. |
|---|
| 14 | # |
|---|
| 15 | # |
|---|
| 16 | # |
|---|
| 17 | |
|---|
| 18 | package provide La 1.0 |
|---|
| 19 | |
|---|
| 20 | |
|---|
| 21 | |
|---|
| 22 | # |
|---|
| 23 | # Package scope and addressed requirements: |
|---|
| 24 | # |
|---|
| 25 | |
|---|
| 26 | # Procedures cover manipulation of vectors and matrices such as |
|---|
| 27 | # scaling, concatenation by rows or columns, subsetting by |
|---|
| 28 | # rows or columns, formatted printing, transpose, dot product, |
|---|
| 29 | # matrix multiplication, solution of linear equation sets, |
|---|
| 30 | # matrix inversion, and eigenvalue/eigenvector solutions. |
|---|
| 31 | # The user can mix vectors and arrays in linear algebra operations. |
|---|
| 32 | # The logic does reasonable conversion of types. |
|---|
| 33 | # Sophisticated operations such as evaluating a custom procedure |
|---|
| 34 | # against each element of a matrix are easily possible. |
|---|
| 35 | # |
|---|
| 36 | # There are typically two calls for a procedure. A plainly |
|---|
| 37 | # named procedure such as "transpose", and another procedure |
|---|
| 38 | # with the suffix "_br" such as "transpose_br". The plainly |
|---|
| 39 | # named procedures expect their data arguments to be passed |
|---|
| 40 | # by value, which is the usual argument passing convention with |
|---|
| 41 | # Tcl programming. The plain calls are designed for ease of |
|---|
| 42 | # interactive use, and in general have been coded to perform |
|---|
| 43 | # more conversion of arguments, more error checking, and to |
|---|
| 44 | # trade off efficiency for convenience. The "_br" procedures |
|---|
| 45 | # are intended for efficent use, with the "_br" indicating that |
|---|
| 46 | # data arguments are passed "By Reference" to avoid copying the |
|---|
| 47 | # data. In Tcl, to pass by reference means that the caller |
|---|
| 48 | # has the data in a named variable, and the caller passes the |
|---|
| 49 | # name of the variable instead of a copy of the data. You can |
|---|
| 50 | # see that passing by reference becomes important for larger |
|---|
| 51 | # vectors and arrays. The _br procedures in general assume |
|---|
| 52 | # that the data arguments have the correct structure, so the |
|---|
| 53 | # caller may need to use the "promote", "demote", and "transpose" |
|---|
| 54 | # calls to prepare arguments for a _br call. |
|---|
| 55 | |
|---|
| 56 | |
|---|
| 57 | |
|---|
| 58 | # |
|---|
| 59 | # La operand formats - scalars, vectors, and arrays |
|---|
| 60 | # |
|---|
| 61 | # Tcl list formats are used for better performance than arrays |
|---|
| 62 | # and better compatibility with C/C++ code. |
|---|
| 63 | # We add dimension information in a simple way at the front of the list. |
|---|
| 64 | # |
|---|
| 65 | # La operands are structured Tcl lists. |
|---|
| 66 | # If the list length is 1 the operand is a scalar, eg., "6.02e23". |
|---|
| 67 | # If the first element in the list is 2, the two following |
|---|
| 68 | # elements provide the number of rows and the number of columns. |
|---|
| 69 | # If the first element is 2: |
|---|
| 70 | # If the number of columns is 0, the operand is a one dimensional |
|---|
| 71 | # vector. So {2 4 0 1 2 3 4} is the vector {1 2 3 4}. Vectors |
|---|
| 72 | # have the same structure as matrices so they can be efficiently |
|---|
| 73 | # modified in place. |
|---|
| 74 | # So a vector of length N has the representation: |
|---|
| 75 | # {2 N 0 v0 v1 v2 ... vN-1} |
|---|
| 76 | # The index into the underlying Tcl list for v[i] is |
|---|
| 77 | # set index [expr 3 + $i] |
|---|
| 78 | # A vector of length N is promoted to an N rows by 1 column |
|---|
| 79 | # matrix if used in an operation where a matrix argument is |
|---|
| 80 | # expected. |
|---|
| 81 | # The transpose of a vector of length N is a 1 row by N columns |
|---|
| 82 | # matrix. |
|---|
| 83 | # If the number of rows and the number of columns are both positive, |
|---|
| 84 | # the operand is a matrix. The data elements follow the number of |
|---|
| 85 | # rows and columns. The data is concatenated as row0, row1, ... |
|---|
| 86 | # so a 2-D matrix, a, has the representation: |
|---|
| 87 | # 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] |
|---|
| 88 | # |
|---|
| 89 | # The index into the underlying Tcl list for a[i,j] is |
|---|
| 90 | # set index [expr 3 + $i * $C + $j] |
|---|
| 91 | # |
|---|
| 92 | # If the first element in the list is > 2, the package assumes the list |
|---|
| 93 | # represents a higher dimensional operand. Logic for higher dimension |
|---|
| 94 | # operands is not currently part of the package. |
|---|
| 95 | # The candidate structure for 3-D data is: |
|---|
| 96 | # 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] |
|---|
| 97 | # P=planes R=rows, C=cols |
|---|
| 98 | # An intuitive view is that the data is multiple 2-D planes of rows |
|---|
| 99 | # and columns listed in order from the 0 plane to the P-1 plane. |
|---|
| 100 | # The index into the underlying Tcl list for a[i,j,k] is |
|---|
| 101 | # set index [expr 4 + $i*$R*$c + $j*$c + $k] |
|---|
| 102 | # |
|---|
| 103 | # |
|---|
| 104 | # |
|---|
| 105 | # La operands, explained by example: |
|---|
| 106 | # |
|---|
| 107 | # <s> a one element scalar (llength == 1) |
|---|
| 108 | # eg., 3.1415 |
|---|
| 109 | # vectors should use the first dimension of 2 |
|---|
| 110 | # 2 1 0 3.1415 ;# a vector of length 1 |
|---|
| 111 | # 2 0 1 3.1415 ;# error - vector should use first dimension only |
|---|
| 112 | # 2 3 0 1 2 3 ;# vector {1 2 3} |
|---|
| 113 | # 2 3 0 1 2 3 4 ;# error - vector has additional element |
|---|
| 114 | # 2 0 0 ;# error - cannot distinguish 1-D or 2-D |
|---|
| 115 | # 2 0 0 <s> ;# error - simpler to have scalars always llength 1 |
|---|
| 116 | # a 2D array of N rows, M columns |
|---|
| 117 | # 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] |
|---|
| 118 | # |
|---|
| 119 | |
|---|
| 120 | # getting started |
|---|
| 121 | # |
|---|
| 122 | # % source "la.tcl" ;# advanced users will use "package require La" |
|---|
| 123 | # % namespace import La::* ;# makes names like "La::show" useable as "show" |
|---|
| 124 | # set x {2 2 3 1 2 3 4 5 6} |
|---|
| 125 | # % show $x |
|---|
| 126 | # 1 2 3 |
|---|
| 127 | # 4 5 6 |
|---|
| 128 | # % show [transpose $x] |
|---|
| 129 | # 1 4 |
|---|
| 130 | # 2 5 |
|---|
| 131 | # 3 6 |
|---|
| 132 | # % set tx [transpose $x] |
|---|
| 133 | # 2 3 2 1 4 2 5 3 6 |
|---|
| 134 | # % show [join_cols $tx $tx] |
|---|
| 135 | # 1 4 1 4 |
|---|
| 136 | # 2 5 2 5 |
|---|
| 137 | # 3 6 3 6 |
|---|
| 138 | # % show [moffset $x 0.2] |
|---|
| 139 | # 1.2 2.2 3.2 |
|---|
| 140 | # 4.2 5.2 6.2 |
|---|
| 141 | # % set v {2 3 0 1 2 3} |
|---|
| 142 | # % show [mmult $x $v] |
|---|
| 143 | # 14.0 |
|---|
| 144 | # 32.0 |
|---|
| 145 | |
|---|
| 146 | |
|---|
| 147 | |
|---|
| 148 | ############################### Package Source Code ##################### |
|---|
| 149 | |
|---|
| 150 | # Advanced Users: |
|---|
| 151 | # The package namespace name is defined here by default. |
|---|
| 152 | # You can specify a different name before using the package. |
|---|
| 153 | global La_namespace |
|---|
| 154 | if { ![info exists La_namespace] } { set La_namespace La } |
|---|
| 155 | |
|---|
| 156 | namespace eval $La_namespace { |
|---|
| 157 | |
|---|
| 158 | namespace export \ |
|---|
| 159 | dim dim_br \ |
|---|
| 160 | dotprod dotprod_br \ |
|---|
| 161 | join_cols join_cols_br \ |
|---|
| 162 | join_rows join_rows_br\ |
|---|
| 163 | lassign lassign_br \ |
|---|
| 164 | madd mdiv mprod msub mat_binary_op mat_binary_op_br \ |
|---|
| 165 | mathprec \ |
|---|
| 166 | mcols mcols_br\ |
|---|
| 167 | mdiag \ |
|---|
| 168 | mdingdong \ |
|---|
| 169 | mevsvd mevsvd_br \ |
|---|
| 170 | mhilbert \ |
|---|
| 171 | mident \ |
|---|
| 172 | mlssvd \ |
|---|
| 173 | mmult mmult_br \ |
|---|
| 174 | mnorms mnorms_br mnormalize mnormalize_br \ |
|---|
| 175 | mrange mrange_br \ |
|---|
| 176 | mround mround_ij \ |
|---|
| 177 | mrows mrows_br\ |
|---|
| 178 | msolve msolve_br \ |
|---|
| 179 | msum msum_br \ |
|---|
| 180 | msvd msvd_br\ |
|---|
| 181 | mat_unary_op madjust moffset mscale mat_unary_op_br \ |
|---|
| 182 | promote promote_br demote demote_br\ |
|---|
| 183 | show show_br \ |
|---|
| 184 | transpose transpose_br \ |
|---|
| 185 | vdiag vdiag_br \ |
|---|
| 186 | vtrim vtrim_br |
|---|
| 187 | |
|---|
| 188 | global La_namespace |
|---|
| 189 | if { $La_namespace != "La" } { |
|---|
| 190 | # help advanced users by moving possible compiled code to |
|---|
| 191 | # their optionally chosen namespace name. |
|---|
| 192 | if { [info commands La::assign_br] == "La::assign_br" } { |
|---|
| 193 | rename La:assign_br $La_namespace::assign_br |
|---|
| 194 | } |
|---|
| 195 | } |
|---|
| 196 | |
|---|
| 197 | ############################# dim #################################### |
|---|
| 198 | ############################# dim #################################### |
|---|
| 199 | ############################# dim #################################### |
|---|
| 200 | # return the dimension of the argument |
|---|
| 201 | # and to some degree verify a proper format |
|---|
| 202 | # {} - null |
|---|
| 203 | # 0 - scalar |
|---|
| 204 | # 1 - vector |
|---|
| 205 | # 2 - array |
|---|
| 206 | # N - higher |
|---|
| 207 | proc dim x { |
|---|
| 208 | return [dim_br x] |
|---|
| 209 | } |
|---|
| 210 | |
|---|
| 211 | proc dim_br name_x { |
|---|
| 212 | upvar $name_x x |
|---|
| 213 | set llen [llength $x] |
|---|
| 214 | if { $llen == 0 } {return {}} |
|---|
| 215 | if { $llen == 1 && [string is double -strict $x] } { return 0 } |
|---|
| 216 | set n [lindex $x 0] |
|---|
| 217 | if { $n == 2 } { |
|---|
| 218 | set d1 [lindex $x 1] |
|---|
| 219 | set d2 [lindex $x 2] |
|---|
| 220 | if {[catch { |
|---|
| 221 | if { $d1 != int($d1) || $d2 != int($d2) } {error {} } |
|---|
| 222 | }] } { error "improper La operand format" } |
|---|
| 223 | if { $d2 == 0 } { |
|---|
| 224 | if { $llen == [expr {$d1 + 3}] } {return 1 } ;# 2 n 0 <n> |
|---|
| 225 | error "improper length of vector" |
|---|
| 226 | } |
|---|
| 227 | if { $llen == [expr {$d1 * $d2 + 3}] } { return 2 } |
|---|
| 228 | error "improper length of matrix" |
|---|
| 229 | } |
|---|
| 230 | if { $n >= 3 && [string is integer -strict $n] } { |
|---|
| 231 | # TBD verify llen as prod of dims + dims + 1 |
|---|
| 232 | return $n |
|---|
| 233 | } |
|---|
| 234 | error "improper La operand format" |
|---|
| 235 | } |
|---|
| 236 | |
|---|
| 237 | |
|---|
| 238 | |
|---|
| 239 | ############################ dotprod ################################# |
|---|
| 240 | ############################ dotprod ################################# |
|---|
| 241 | ############################ dotprod ################################# |
|---|
| 242 | # |
|---|
| 243 | # Dot Product = sum over i, Ai * Bi |
|---|
| 244 | # |
|---|
| 245 | # Can work columns or rows in matrices |
|---|
| 246 | # because indexing increment is passed. |
|---|
| 247 | # Default arguments work with vectors, Nx1, or 1XN matrices |
|---|
| 248 | # |
|---|
| 249 | proc dotprod {a b {N {}} {a0index 3} {b0index 3} {ainc 1} {binc 1}} { |
|---|
| 250 | return [dotprod_br a b $N $a0index $b0index $ainc $binc] |
|---|
| 251 | } |
|---|
| 252 | |
|---|
| 253 | # |
|---|
| 254 | # always returns a scalar result |
|---|
| 255 | # |
|---|
| 256 | proc dotprod_br {a_in b_in {N {}} {a0index 3} {b0index 3} {ainc 1} {binc 1}} { |
|---|
| 257 | upvar $a_in a |
|---|
| 258 | upvar $b_in b |
|---|
| 259 | set z 0.0 |
|---|
| 260 | if { $N == {}} { ;# default length assumes vector, 1xN, or Nx1 first argument |
|---|
| 261 | set len [llength $a] |
|---|
| 262 | set N [expr {$len-3}] |
|---|
| 263 | # and then see if that matches one of b's dimensions |
|---|
| 264 | if { [lindex $b 1] != $N && [lindex $b 2] != $N } { |
|---|
| 265 | error "Assumed length does not seem proper." |
|---|
| 266 | } |
|---|
| 267 | } |
|---|
| 268 | for {set i 0} {$i < $N} {incr i} { |
|---|
| 269 | set aval [lindex $a $a0index] |
|---|
| 270 | set bval [lindex $b $b0index] |
|---|
| 271 | set z [expr { $z + $aval * $bval}] |
|---|
| 272 | incr a0index $ainc |
|---|
| 273 | incr b0index $binc |
|---|
| 274 | } |
|---|
| 275 | return $z |
|---|
| 276 | } |
|---|
| 277 | |
|---|
| 278 | |
|---|
| 279 | ########################### join_cols ################################ |
|---|
| 280 | ########################### join_cols ################################ |
|---|
| 281 | ########################### join_cols ################################ |
|---|
| 282 | # |
|---|
| 283 | # combine vector or matrix args as added columns |
|---|
| 284 | # |
|---|
| 285 | proc join_cols {a b} { |
|---|
| 286 | set dimA [dim_br a] |
|---|
| 287 | set dimB [dim_br b] |
|---|
| 288 | if { $dimA > 2 || $dimB > 2 } { |
|---|
| 289 | error "cannot handle > 2D data" |
|---|
| 290 | } |
|---|
| 291 | if { $dimA < 2 } {promote_br a} |
|---|
| 292 | if { $dimB < 2 } {promote_br b} |
|---|
| 293 | join_cols_br a b |
|---|
| 294 | return $a |
|---|
| 295 | } |
|---|
| 296 | |
|---|
| 297 | |
|---|
| 298 | proc join_cols_br {a_in b_in {c_out {}}} { |
|---|
| 299 | upvar $a_in a |
|---|
| 300 | upvar $b_in b |
|---|
| 301 | if { $c_out == {}} {set c_out $a_in } |
|---|
| 302 | upvar $c_out c |
|---|
| 303 | # have to match row counts |
|---|
| 304 | set rows [lindex $a 1] |
|---|
| 305 | set acols [lindex $a 2] |
|---|
| 306 | set bcols [lindex $b 2] |
|---|
| 307 | if { [lindex $b 1] != $rows } { |
|---|
| 308 | error "cannot append columns with inequal rows A\[$rows,\] + B\[[lindex $b 1],\]" |
|---|
| 309 | } |
|---|
| 310 | set total_cols [expr {$acols + $bcols}] |
|---|
| 311 | set result [list 2 $rows $total_cols] |
|---|
| 312 | for {set row 0} {$row < $rows} {incr row} { |
|---|
| 313 | for {set col 0} {$col < $acols} {incr col} { |
|---|
| 314 | set i_a [expr {3 + $row * $acols + $col}] |
|---|
| 315 | lappend result [lindex $a $i_a] |
|---|
| 316 | } |
|---|
| 317 | for {set col 0} {$col < $bcols} {incr col} { |
|---|
| 318 | set i_b [expr {3 + $row * $bcols + $col}] |
|---|
| 319 | lappend result [lindex $b $i_b] |
|---|
| 320 | } |
|---|
| 321 | } |
|---|
| 322 | set c $result |
|---|
| 323 | return $c_out |
|---|
| 324 | } |
|---|
| 325 | |
|---|
| 326 | ########################### join_rows ################################ |
|---|
| 327 | ########################### join_rows ################################ |
|---|
| 328 | ########################### join_rows ################################ |
|---|
| 329 | # |
|---|
| 330 | # combine vector or matrix args as added rows |
|---|
| 331 | # |
|---|
| 332 | proc join_rows {a b} { |
|---|
| 333 | set dimA [dim_br a] |
|---|
| 334 | set dimB [dim_br b] |
|---|
| 335 | if { $dimA > 2 || $dimB > 2 } { |
|---|
| 336 | error "cannot handle > 2D data" |
|---|
| 337 | } |
|---|
| 338 | if { $dimA < 2 } {promote_br a} |
|---|
| 339 | if { $dimB < 2 } {promote_br b} |
|---|
| 340 | join_rows_br a b |
|---|
| 341 | return $a |
|---|
| 342 | } |
|---|
| 343 | |
|---|
| 344 | |
|---|
| 345 | proc join_rows_br {a_in b_in {c_out {}}} { |
|---|
| 346 | upvar $a_in a |
|---|
| 347 | upvar $b_in b |
|---|
| 348 | if { $c_out == {}} {set c_out $a_in } |
|---|
| 349 | upvar $c_out c |
|---|
| 350 | # have to match col counts |
|---|
| 351 | set arows [lindex $a 1] |
|---|
| 352 | set brows [lindex $b 1] |
|---|
| 353 | set acols [lindex $a 2] |
|---|
| 354 | set bcols [lindex $b 2] |
|---|
| 355 | if { $acols != $bcols } { |
|---|
| 356 | error "cannot append rows with inequal columns A\[,$acols\] + B\[,$bcols\]" |
|---|
| 357 | } |
|---|
| 358 | set total_rows [expr {$arows + $brows}] |
|---|
| 359 | set c [concat 2 $total_rows $acols [lrange $a 3 end] [lrange $b 3 end]] |
|---|
| 360 | return $c_out |
|---|
| 361 | } |
|---|
| 362 | |
|---|
| 363 | ########################### lassign ################################## |
|---|
| 364 | ########################### lassign ################################## |
|---|
| 365 | ########################### lassign ################################## |
|---|
| 366 | |
|---|
| 367 | # |
|---|
| 368 | # replace a single value in a Tcl list |
|---|
| 369 | # This call assumes the caller knows the desired position in the list. |
|---|
| 370 | # |
|---|
| 371 | proc lassign {x index value} { |
|---|
| 372 | return [lreplace $x $index $index $value] |
|---|
| 373 | } |
|---|
| 374 | |
|---|
| 375 | # the mysterious K procedure improves efficiency because of how list |
|---|
| 376 | # manipulation is performed (thank you Kevin Kenny) |
|---|
| 377 | proc K {x y} {set x} |
|---|
| 378 | |
|---|
| 379 | # prototype list element replace by value - this can be done in C code |
|---|
| 380 | # without copying the list so it can be efficient, especially since Tcl |
|---|
| 381 | # will maintain an internal representation of a list already |
|---|
| 382 | # split into elements. |
|---|
| 383 | # lreplace forces you to copy the whole list. |
|---|
| 384 | if { [info commands ${La_namespace}::lassign_br] == {}} { |
|---|
| 385 | proc lassign_br {listname index value} { |
|---|
| 386 | upvar $listname thelist |
|---|
| 387 | set thelist [lreplace [K $thelist [set thelist {}]] $index $index $value] |
|---|
| 388 | return $listname |
|---|
| 389 | } |
|---|
| 390 | } |
|---|
| 391 | |
|---|
| 392 | |
|---|
| 393 | ################ madd,mdiv, mprod, msub, mat_binary_op ############## |
|---|
| 394 | ################ madd,mdiv, mprod, msub, mat_binary_op ############## |
|---|
| 395 | ################ madd,mdiv, mprod, msub, mat_binary_op ############## |
|---|
| 396 | # |
|---|
| 397 | # vector/matrix binary operations on elements like addition |
|---|
| 398 | # |
|---|
| 399 | proc mat_binary_op {a b op} { |
|---|
| 400 | set dimA [dim_br a] |
|---|
| 401 | set dimB [dim_br b] |
|---|
| 402 | if { $dimA > 2 || $dimB > 2 } { |
|---|
| 403 | error "cannot handle > 2D data" |
|---|
| 404 | } |
|---|
| 405 | if { $dimA < 1 || $dimA < $dimB } {promote_br a} |
|---|
| 406 | if { $dimB < 1 || $dimB < $dimA } {promote_br b} |
|---|
| 407 | mat_binary_op_br a b $op c |
|---|
| 408 | return $c |
|---|
| 409 | } |
|---|
| 410 | |
|---|
| 411 | # matrix addition and similar |
|---|
| 412 | proc madd {a b} { return [mat_binary_op $a $b +] } |
|---|
| 413 | proc mdiv {a b} { return [mat_binary_op $a $b {*1.0/}] } |
|---|
| 414 | proc mprod {a b} { return [mat_binary_op $a $b *] } |
|---|
| 415 | proc msub {a b} { return [mat_binary_op $a $b -] } |
|---|
| 416 | |
|---|
| 417 | |
|---|
| 418 | # perform a binary operation using corresponding elements of |
|---|
| 419 | # two matrices or vectors - op is commonly +, -, *, / |
|---|
| 420 | # but you can also add formula constructions as in |
|---|
| 421 | # " *2.0 + 3.4*" |
|---|
| 422 | # Be careful with division since Tcl tries to perform integer |
|---|
| 423 | # division. You may want to use "*1.0 / " to insure floating point math |
|---|
| 424 | # |
|---|
| 425 | proc mat_binary_op_br {a_in b_in op {c_out {}}} { |
|---|
| 426 | upvar $a_in a |
|---|
| 427 | upvar $b_in b |
|---|
| 428 | if { $c_out == {}} { set c_out $a_in } |
|---|
| 429 | upvar $c_out c |
|---|
| 430 | set NR [lindex $a 1] |
|---|
| 431 | set MC [lindex $a 2] |
|---|
| 432 | set NR2 [lindex $b 1] |
|---|
| 433 | set MC2 [lindex $b 2] |
|---|
| 434 | if { $NR != $NR2 || $MC != $MC2 } { |
|---|
| 435 | error "arguments are not conformable A\[$NR,$MC\] vs B\[$NR2,$MC2\]" |
|---|
| 436 | } |
|---|
| 437 | set result [list 2 $NR $MC] |
|---|
| 438 | set imax [llength $a] |
|---|
| 439 | for {set i 3} {$i < $imax} {incr i} { |
|---|
| 440 | set ai [lindex $a $i] |
|---|
| 441 | set bi [lindex $b $i] |
|---|
| 442 | lappend result [expr $ai $op $bi] |
|---|
| 443 | } |
|---|
| 444 | set c $result |
|---|
| 445 | return $c_out |
|---|
| 446 | } |
|---|
| 447 | |
|---|
| 448 | |
|---|
| 449 | ########################## mathprec ################################ |
|---|
| 450 | ########################## mathprec ################################ |
|---|
| 451 | ########################## mathprec ################################ |
|---|
| 452 | # |
|---|
| 453 | # test math precision |
|---|
| 454 | |
|---|
| 455 | # what is the smallest number such that |
|---|
| 456 | # 1+epsilon > 1 |
|---|
| 457 | # does the machine truncate or round-off |
|---|
| 458 | # what is the radix and number of digits |
|---|
| 459 | # |
|---|
| 460 | # Intel PC (IEEE 8 byte floating point) results: |
|---|
| 461 | # radix=2.0 digits=53 epsilon=2.22044604925e-016 method=truncation |
|---|
| 462 | # |
|---|
| 463 | proc mathprec {{puts puts}} { |
|---|
| 464 | for {set test 1.0} {$test + 1.0 != $test} {set test [expr {$test*2.0}]} { } |
|---|
| 465 | for {set diff 1.0} {$test + $diff == $test} {set diff [expr {$diff + 1.0}]} { } |
|---|
| 466 | set radix [expr {($test + $diff) - $test}] |
|---|
| 467 | if { $diff < $radix } { set method rounding } else {set method truncation } |
|---|
| 468 | set digits 0 |
|---|
| 469 | for {set test 1.0} {$test + 1 != $test} {set test [expr {$test*$radix}]} { |
|---|
| 470 | incr digits |
|---|
| 471 | } |
|---|
| 472 | set epsilon [expr {pow($radix,(1.0-$digits))}] |
|---|
| 473 | if {$puts != {}} { |
|---|
| 474 | $puts "radix=$radix digits=$digits epsilon=$epsilon method=$method" |
|---|
| 475 | } |
|---|
| 476 | return $epsilon |
|---|
| 477 | } |
|---|
| 478 | |
|---|
| 479 | ############################ mcols ################################# |
|---|
| 480 | ############################ mcols ################################# |
|---|
| 481 | ############################ mcols ################################# |
|---|
| 482 | # |
|---|
| 483 | # you can use mcols and mrows to determine the number of columns or |
|---|
| 484 | # rows and keep your code isolated from the details of the |
|---|
| 485 | # data representation |
|---|
| 486 | # |
|---|
| 487 | proc mcols m { return [lindex $m 2] } |
|---|
| 488 | proc mcols_br name_m { upvar $name_m m ; return [lindex $m 2] } |
|---|
| 489 | proc mrows m { return [lindex $m 1] } |
|---|
| 490 | proc mrows_br name_m { upvar $name_m m ; return [lindex $m 1] } |
|---|
| 491 | |
|---|
| 492 | ############################ mdiag ################################# |
|---|
| 493 | ############################ mdiag ################################# |
|---|
| 494 | ############################ mdiag ################################# |
|---|
| 495 | # |
|---|
| 496 | # create a diagonal matrix from a vector, an Nx1 matrix, or a 1XN matrix |
|---|
| 497 | # |
|---|
| 498 | proc mdiag v { |
|---|
| 499 | set dim [dim_br v] |
|---|
| 500 | if { $dim == 2 } { demote_br v } |
|---|
| 501 | if { [lindex $v 2] != 0 } { error "expected vector argument" } |
|---|
| 502 | set N [lindex $v 1] |
|---|
| 503 | set result [list 2 $N $N] |
|---|
| 504 | for {set row 0} {$row < $N} {incr row} { |
|---|
| 505 | for {set col 0} {$col < $N} {incr col} { |
|---|
| 506 | set iv [expr {3 + $row}] |
|---|
| 507 | if { $row == $col } { lappend result [lindex $v $iv] }\ |
|---|
| 508 | else { lappend result 0 } |
|---|
| 509 | } |
|---|
| 510 | } |
|---|
| 511 | return $result |
|---|
| 512 | } |
|---|
| 513 | |
|---|
| 514 | ############################ mdingdong ############################## |
|---|
| 515 | ############################ mdingdong ############################## |
|---|
| 516 | ############################ mdingdong ############################## |
|---|
| 517 | # |
|---|
| 518 | # create the Ding Dong test matrix, a Cauchy matrix that is |
|---|
| 519 | # represented inexactly in the machine, but very stable for |
|---|
| 520 | # inversion by elimination methods |
|---|
| 521 | # |
|---|
| 522 | proc mdingdong N { |
|---|
| 523 | if { ![string is integer -strict $N] } { error "improper size \"$N\"" } |
|---|
| 524 | set result [list 2 $N $N] |
|---|
| 525 | for {set row 0} {$row < $N} {incr row} { |
|---|
| 526 | for {set col 0} {$col < $N} {incr col} { |
|---|
| 527 | lappend result [expr {0.5/($N - $col - $row - 0.5)}] |
|---|
| 528 | } |
|---|
| 529 | } |
|---|
| 530 | return $result |
|---|
| 531 | } |
|---|
| 532 | |
|---|
| 533 | ############################ mevsvd ################################# |
|---|
| 534 | ############################ mevsvd ################################# |
|---|
| 535 | ############################ mevsvd ################################# |
|---|
| 536 | # |
|---|
| 537 | # eigenvectors/eigenvalues of a real symmetric matrix by |
|---|
| 538 | # singular value decomposition |
|---|
| 539 | # |
|---|
| 540 | proc mevsvd {A {eps 2.3e-16}} { |
|---|
| 541 | mevsvd_br A evals |
|---|
| 542 | puts "eigenvalues=[show $evals]" |
|---|
| 543 | return $A |
|---|
| 544 | } |
|---|
| 545 | |
|---|
| 546 | |
|---|
| 547 | # the eigenvectors are returned in place as the columns of A |
|---|
| 548 | # the eigenvalues are returned as a vector |
|---|
| 549 | |
|---|
| 550 | proc mevsvd_br {A_in_out evals_out {eps 2.3e-16}} { |
|---|
| 551 | upvar $A_in_out A |
|---|
| 552 | upvar $evals_out evals |
|---|
| 553 | set n [lindex $A 1] |
|---|
| 554 | if { [lindex $A 2] != $n } { error "expecting square matrix" } |
|---|
| 555 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 556 | set ii [expr {3 + $i*$n + $i}] |
|---|
| 557 | set v [lindex $A $ii] |
|---|
| 558 | for {set j 0} {$j < $n} {incr j} { |
|---|
| 559 | if { $i != $j } { |
|---|
| 560 | set ij [expr {3 + $i*$n + $j}] |
|---|
| 561 | set Aij [lindex $A $ij] |
|---|
| 562 | set v [expr {$v - abs($Aij)}] |
|---|
| 563 | } |
|---|
| 564 | } |
|---|
| 565 | if { ![info exists h] } { set h $v }\ |
|---|
| 566 | elseif { $v < $h } { set h $v } |
|---|
| 567 | } |
|---|
| 568 | # h is lower bound on Gershgorin region |
|---|
| 569 | if { $h <= $eps } { |
|---|
| 570 | set h [expr {$h - sqrt($eps)}] |
|---|
| 571 | # try to make smallest eigenvalue positive and not too small |
|---|
| 572 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 573 | set ii [expr {3 + $i*$n + $i}] |
|---|
| 574 | set Aii [lindex $A $ii] |
|---|
| 575 | lassign_br A $ii [expr {$Aii - $h}] |
|---|
| 576 | } |
|---|
| 577 | }\ |
|---|
| 578 | else { |
|---|
| 579 | set h 0.0 |
|---|
| 580 | } |
|---|
| 581 | |
|---|
| 582 | set count 0 |
|---|
| 583 | # top of the iteration |
|---|
| 584 | for {set isweep 0} {$isweep < 30 && $count < $n*($n-1)/2} {incr isweep} { |
|---|
| 585 | set count 0 ;# count of rotations in a sweep |
|---|
| 586 | |
|---|
| 587 | for {set j 0} {$j < [expr {$n-1}]} {incr j} { |
|---|
| 588 | for {set k [expr {$j+1}]} {$k < $n} {incr k} { |
|---|
| 589 | set p [set q [set r 0.0]] |
|---|
| 590 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 591 | set ij [expr {3+$i*$n+$j}] |
|---|
| 592 | set ik [expr {3+$i*$n+$k}] |
|---|
| 593 | set Aij [lindex $A $ij] |
|---|
| 594 | set Aik [lindex $A $ik] |
|---|
| 595 | set p [expr {$p + $Aij*$Aik}] |
|---|
| 596 | set q [expr {$q + $Aij*$Aij}] |
|---|
| 597 | set r [expr {$r + $Aik*$Aik}] |
|---|
| 598 | } |
|---|
| 599 | if { 1.0 >= 1.0 + abs($p/sqrt($q*$r)) } { |
|---|
| 600 | if { $q >= $r } { |
|---|
| 601 | incr count |
|---|
| 602 | # no rotation needed |
|---|
| 603 | continue |
|---|
| 604 | } |
|---|
| 605 | } |
|---|
| 606 | set q [expr {$q-$r}] |
|---|
| 607 | set v [expr {sqrt(4.0*$p*$p + $q*$q)}] |
|---|
| 608 | if { $v == 0.0 } continue |
|---|
| 609 | if { $q >= 0.0 } { |
|---|
| 610 | set c [expr {sqrt(($v+$q)/(2.0*$v))}] |
|---|
| 611 | set s [expr {$p/($v*$c)}] |
|---|
| 612 | }\ |
|---|
| 613 | else { |
|---|
| 614 | set s [expr {sqrt(($v-$q)/(2.0*$v))}] |
|---|
| 615 | if { $p < 0.0 } { set s [expr {0.0-$s}] } |
|---|
| 616 | set c [expr {$p/($v*$s)}] |
|---|
| 617 | } |
|---|
| 618 | # rotation |
|---|
| 619 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 620 | set ij [expr {3+$i*$n+$j}] |
|---|
| 621 | set ik [expr {3+$i*$n+$k}] |
|---|
| 622 | set Aij [lindex $A $ij] |
|---|
| 623 | set Aik [lindex $A $ik] |
|---|
| 624 | lassign_br A $ij [expr {$Aij*$c + $Aik*$s}] |
|---|
| 625 | lassign_br A $ik [expr {$Aik*$c - $Aij*$s}] |
|---|
| 626 | } |
|---|
| 627 | } ;# k |
|---|
| 628 | } ;# j |
|---|
| 629 | #puts "pass=$isweep skipped rotations=$count" |
|---|
| 630 | } ;# isweep |
|---|
| 631 | |
|---|
| 632 | # now columns are orthogonal, rescale |
|---|
| 633 | # and flip signs if all negative or zero |
|---|
| 634 | set evals [list 2 $n 0] |
|---|
| 635 | for {set j 0} {$j < $n} {incr j} { |
|---|
| 636 | set s 0.0 |
|---|
| 637 | set notpositive 0 |
|---|
| 638 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 639 | set ij [expr {3+$i*$n+$j}] |
|---|
| 640 | set Aij [lindex $A $ij] |
|---|
| 641 | if { $Aij <= 0.0 } { incr notpositive } |
|---|
| 642 | set s [expr {$s + $Aij*$Aij}] |
|---|
| 643 | } |
|---|
| 644 | set s [expr {sqrt($s)}] |
|---|
| 645 | if { $notpositive == $n } { set sf [expr {0.0-$s}] }\ |
|---|
| 646 | else { set sf $s } |
|---|
| 647 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 648 | set ij [expr {3+$i*$n+$j}] |
|---|
| 649 | set Aij [lindex $A $ij] |
|---|
| 650 | lassign_br A $ij [expr {$Aij/$sf}] |
|---|
| 651 | } |
|---|
| 652 | lappend evals [expr {$s+$h}] |
|---|
| 653 | } |
|---|
| 654 | return $A_in_out |
|---|
| 655 | } |
|---|
| 656 | |
|---|
| 657 | ############################ mhilbert ############################### |
|---|
| 658 | ############################ mhilbert ############################### |
|---|
| 659 | ############################ mhilbert ############################### |
|---|
| 660 | # |
|---|
| 661 | # create the Hilbert test matrix which is notorious for being |
|---|
| 662 | # ill conditioned for eigenvector/eigenvalue solutions |
|---|
| 663 | # |
|---|
| 664 | proc mhilbert N { |
|---|
| 665 | if { ![string is integer -strict $N] } { error "improper size \"$N\"" } |
|---|
| 666 | set result [list 2 $N $N] |
|---|
| 667 | for {set row 0} {$row < $N} {incr row} { |
|---|
| 668 | for {set col 0} {$col < $N} {incr col} { |
|---|
| 669 | lappend result [expr {1.0/($col + $row +1.0)}] |
|---|
| 670 | } |
|---|
| 671 | } |
|---|
| 672 | return $result |
|---|
| 673 | } |
|---|
| 674 | |
|---|
| 675 | ############################ mident ################################# |
|---|
| 676 | ############################ mident ################################# |
|---|
| 677 | ############################ mident ################################# |
|---|
| 678 | # |
|---|
| 679 | # create an identity matrix of order N |
|---|
| 680 | # |
|---|
| 681 | proc mident N { |
|---|
| 682 | if { ![string is integer -strict $N] } { error "improper size \"$N\"" } |
|---|
| 683 | set result [list 2 $N $N] |
|---|
| 684 | for {set row 0} {$row < $N} {incr row} { |
|---|
| 685 | for {set col 0} {$col < $N} {incr col} { |
|---|
| 686 | if { $row == $col } { lappend result 1 }\ |
|---|
| 687 | else { lappend result 0 } |
|---|
| 688 | } |
|---|
| 689 | } |
|---|
| 690 | return $result |
|---|
| 691 | } |
|---|
| 692 | |
|---|
| 693 | ############################### mlssvd ######################### |
|---|
| 694 | ############################### mlssvd ######################### |
|---|
| 695 | ############################### mlssvd ######################### |
|---|
| 696 | # |
|---|
| 697 | # solve the over-determined linear equations in the least squares |
|---|
| 698 | # sense using SVD |
|---|
| 699 | # |
|---|
| 700 | # Ax ~ y |
|---|
| 701 | # y[m] is dependent variable such as a measured outcome |
|---|
| 702 | # x[n] vector of independent variables |
|---|
| 703 | # A[m,n] - each row is a set dependent variable values, |
|---|
| 704 | # the first column is usually 1 to allow for a constant |
|---|
| 705 | # in the regression |
|---|
| 706 | # |
|---|
| 707 | # q is the minimum singular value, lessor values are treated as 0 |
|---|
| 708 | # |
|---|
| 709 | proc mlssvd {A y {q 0.0} {puts puts} {epsilon 2.3e-16}} { |
|---|
| 710 | # A[m,n] |
|---|
| 711 | set m [lindex $A 1] |
|---|
| 712 | set n [lindex $A 2] |
|---|
| 713 | promote_br y ;# now expect y[m,1] |
|---|
| 714 | if { [lindex $y 1] != $m } { |
|---|
| 715 | error "cannot conform A\[$m,$n\]*X\[$n] = y\[[lindex $y 1],[lindex $y 2]]" |
|---|
| 716 | } |
|---|
| 717 | msvd_br A S V |
|---|
| 718 | # now A has been transformed to U[m,n] |
|---|
| 719 | # S[n], V[n,n] |
|---|
| 720 | if { $puts != {}} { |
|---|
| 721 | $puts singular\ values=[show $S %.6g] |
|---|
| 722 | } |
|---|
| 723 | set tol [expr {$epsilon * $epsilon * $n * $n}] |
|---|
| 724 | # form Utrans*y into g |
|---|
| 725 | set g [list 2 $n 0] |
|---|
| 726 | for {set j 0} {$j < $n} {incr j} { |
|---|
| 727 | set s 0.0 |
|---|
| 728 | for {set i 0} {$i < $m} {incr i} { |
|---|
| 729 | set ij [expr {3 + $i*$n + $j}] |
|---|
| 730 | set Aij [lindex $A $ij] |
|---|
| 731 | set yi [lindex $y [expr {3 + $i}]] |
|---|
| 732 | set s [expr {$s + $Aij*$yi}] |
|---|
| 733 | } |
|---|
| 734 | lappend g $s ;# g[j] = $s |
|---|
| 735 | } |
|---|
| 736 | # form VS+g = VS+Utrans*g |
|---|
| 737 | set x [list 2 $n 0] |
|---|
| 738 | for {set j 0} {$j < $n} {incr j} { |
|---|
| 739 | set s 0.0 |
|---|
| 740 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 741 | set iindex [expr {$i+3}] |
|---|
| 742 | set zi [lindex $S $iindex] |
|---|
| 743 | if { $zi > $q } { |
|---|
| 744 | set ji [expr {3 + $j*$n+$i}] |
|---|
| 745 | set Vji [lindex $V $ji] |
|---|
| 746 | set gi [lindex $g $iindex] |
|---|
| 747 | set s [expr {$s + $Vji*$gi/$zi}] |
|---|
| 748 | } |
|---|
| 749 | } |
|---|
| 750 | lappend x $s |
|---|
| 751 | } |
|---|
| 752 | return $x |
|---|
| 753 | } |
|---|
| 754 | |
|---|
| 755 | ############################ mmult ################################## |
|---|
| 756 | ############################ mmult ################################## |
|---|
| 757 | ############################ mmult ################################## |
|---|
| 758 | # |
|---|
| 759 | # matrix multiplication |
|---|
| 760 | # |
|---|
| 761 | # A[p,q] x B[q,r] = C[p,r] |
|---|
| 762 | # Cij = dot product A[i,] row x B[,j] col |
|---|
| 763 | # Cij = sum Aik * Bkj ; k=0..q-1 |
|---|
| 764 | # |
|---|
| 765 | # rules: Vector arguments are always promoted to Nx1 arrays. |
|---|
| 766 | # So chances are if you are using one as a left operand |
|---|
| 767 | # you probably intend to use the transpose of it (1xN), |
|---|
| 768 | # which is easily done using the transpose procedure. |
|---|
| 769 | # |
|---|
| 770 | proc mmult {a b} { |
|---|
| 771 | set dimA [dim_br a] |
|---|
| 772 | set dimB [dim_br b] |
|---|
| 773 | if { $dimA > 2 || $dimB > 2 } { |
|---|
| 774 | error "cannot handle > 2D data" |
|---|
| 775 | } |
|---|
| 776 | if { $dimA < 2 } {promote_br a} |
|---|
| 777 | if { $dimB < 2 } {promote_br b} |
|---|
| 778 | mmult_br a b c |
|---|
| 779 | return $c |
|---|
| 780 | } |
|---|
| 781 | |
|---|
| 782 | # caller is responsible to provide 2D arguments |
|---|
| 783 | proc mmult_br {a_in b_in c_out} { |
|---|
| 784 | upvar $a_in a |
|---|
| 785 | upvar $b_in b |
|---|
| 786 | upvar $c_out c |
|---|
| 787 | set p [lindex $a 1] |
|---|
| 788 | set q [lindex $a 2] |
|---|
| 789 | set q2 [lindex $b 1] |
|---|
| 790 | set r [lindex $b 2] |
|---|
| 791 | if { $q2 != $q } { |
|---|
| 792 | error "matrices are not conformable A\[$p,$q\] x B\[$q2,$r\]" |
|---|
| 793 | } |
|---|
| 794 | # Cij = dot product A[i,] row x B[,j] col |
|---|
| 795 | set c [list 2 $p $r] |
|---|
| 796 | set a_base 3 |
|---|
| 797 | set ainc 1 |
|---|
| 798 | set binc $r |
|---|
| 799 | for {set row 0} {$row < $p} {incr row} { |
|---|
| 800 | set b_base 3 |
|---|
| 801 | for {set col 0} {$col < $r} {incr col} { |
|---|
| 802 | lappend c [dotprod_br a b $q $a_base $b_base $ainc $binc] |
|---|
| 803 | incr b_base 1 |
|---|
| 804 | } |
|---|
| 805 | incr a_base $q |
|---|
| 806 | } |
|---|
| 807 | |
|---|
| 808 | return $c_out |
|---|
| 809 | } |
|---|
| 810 | |
|---|
| 811 | ######################## mnorms, mnormalize ########################## |
|---|
| 812 | ######################## mnorms, mnormalize ########################## |
|---|
| 813 | ######################## mnorms, mnormalize ########################## |
|---|
| 814 | # |
|---|
| 815 | # compute the means and sigmas of each column |
|---|
| 816 | # |
|---|
| 817 | # |
|---|
| 818 | proc mnorms a { |
|---|
| 819 | mnorms_br a means sigmas |
|---|
| 820 | return [transpose [join_cols $means $sigmas]] |
|---|
| 821 | } |
|---|
| 822 | |
|---|
| 823 | proc mnorms_br {a_in means_out sigmas_out} { |
|---|
| 824 | upvar $a_in a |
|---|
| 825 | if { [dim_br a] != 2 } { error "expecting matrix" } |
|---|
| 826 | set NR [lindex $a 1] |
|---|
| 827 | if { $NR < 2 } { error "insufficient rows to calculate sigmas" } |
|---|
| 828 | set NC [lindex $a 2] |
|---|
| 829 | # vector results are returned |
|---|
| 830 | set means [set sigmas [list 2 $NC 0]] |
|---|
| 831 | for {set i 0} {$i < $NC} {incr i} { |
|---|
| 832 | mrange_br a colvect $i $i |
|---|
| 833 | # mean, stddev are Hume C-code Tcl commands |
|---|
| 834 | set coldata [lrange $colvect 3 end] |
|---|
| 835 | lappend means [mean $coldata] |
|---|
| 836 | lappend sigmas [stddev $coldata] |
|---|
| 837 | } |
|---|
| 838 | upvar $means_out m |
|---|
| 839 | set m $means |
|---|
| 840 | upvar $sigmas_out s |
|---|
| 841 | set s $sigmas |
|---|
| 842 | return [list $means_out $sigmas_out] |
|---|
| 843 | } |
|---|
| 844 | |
|---|
| 845 | # |
|---|
| 846 | # normalize each column by subtracting the corresponding mean and then |
|---|
| 847 | # dividing by the corresponding sigma |
|---|
| 848 | # |
|---|
| 849 | |
|---|
| 850 | proc mnormalize {a means sigmas} { |
|---|
| 851 | if { [dim_br a] == 1 } { promote_br a } |
|---|
| 852 | mnormalize_br a means sigmas |
|---|
| 853 | return $a |
|---|
| 854 | } |
|---|
| 855 | |
|---|
| 856 | # |
|---|
| 857 | # the code will work with means and sigmas as vectors, Nx1, or 1xN matrices |
|---|
| 858 | # |
|---|
| 859 | proc mnormalize_br {a_in means_in sigmas_in {c_out {}}} { |
|---|
| 860 | upvar $a_in a |
|---|
| 861 | upvar $means_in means |
|---|
| 862 | upvar $sigmas_in sigmas |
|---|
| 863 | if { [dim_br a] != 2 } { error "expecting matrix" } |
|---|
| 864 | set NR [lindex $a 1] |
|---|
| 865 | set NC [lindex $a 2] |
|---|
| 866 | set nm [llength $means] |
|---|
| 867 | set ns [llength $sigmas] |
|---|
| 868 | if { $nm != $ns || $nm != [expr 3+$NC] } { |
|---|
| 869 | error "non-conformable arguments a[$NR,$NC] means[expr {$nm-3}] sigmas[expr {$ns-3}]" |
|---|
| 870 | } |
|---|
| 871 | set result [list 2 $NR $NC] |
|---|
| 872 | for {set row 0} {$row < $NR} {incr row} { |
|---|
| 873 | for {set i 0} {$i < $NC} {incr i} { |
|---|
| 874 | set vindex [expr {3 + $i}] |
|---|
| 875 | set mean [lindex $means $vindex] |
|---|
| 876 | set sigma [lindex $sigmas $vindex] |
|---|
| 877 | set aindex [expr {3 + $row*$NC + $i}] |
|---|
| 878 | set val [lindex $a $aindex] |
|---|
| 879 | lappend result [expr { (0.0 + $val - $mean)/$sigma }] |
|---|
| 880 | } |
|---|
| 881 | } |
|---|
| 882 | if { $c_out == {}} { set c_out $a_in } |
|---|
| 883 | upvar $c_out c |
|---|
| 884 | set c $result |
|---|
| 885 | return $c_out |
|---|
| 886 | } |
|---|
| 887 | |
|---|
| 888 | # the Hume dmh_wish has these in C code, here they are for other shells |
|---|
| 889 | if { [info commands mean] != "mean" } { |
|---|
| 890 | proc ::mean numlist { |
|---|
| 891 | set len [llength $numlist] |
|---|
| 892 | if { $len == 0 } { return {}} |
|---|
| 893 | if { $len == 1 } { return [expr {double($numlist)}] } |
|---|
| 894 | set s 0.0 |
|---|
| 895 | for {set i 0} {$i < $len} {incr i} { |
|---|
| 896 | set x [lindex $numlist $i] |
|---|
| 897 | set s [expr {$s + $x}] |
|---|
| 898 | } |
|---|
| 899 | return [expr {$s/$len}] |
|---|
| 900 | } |
|---|
| 901 | } |
|---|
| 902 | if { [info commands stddev] != "stddev" } { |
|---|
| 903 | proc ::stddev numlist { |
|---|
| 904 | set len [llength $numlist] |
|---|
| 905 | if { $len < 2 } { return {} } |
|---|
| 906 | set sx [set sxx 0.0] |
|---|
| 907 | for {set i 0} {$i < $len} {incr i} { |
|---|
| 908 | set x [lindex $numlist $i] |
|---|
| 909 | set sx [expr {$sx + $x}] |
|---|
| 910 | set sxx [expr {$sxx + $x * $x}] |
|---|
| 911 | } |
|---|
| 912 | # the real world is full of surprises like stdev {2.41 2.41 2.41} |
|---|
| 913 | # causing a domain error, so you get strange code |
|---|
| 914 | set diff [expr {$sxx - $sx*$sx/$len}] |
|---|
| 915 | if { $diff < 0.0 || ($diff + 0.0 != $diff) } {return 0.0 } |
|---|
| 916 | return [expr {sqrt($diff/($len-1.0))}] |
|---|
| 917 | } |
|---|
| 918 | } |
|---|
| 919 | |
|---|
| 920 | #################### mrange ########################################## |
|---|
| 921 | #################### mrange ########################################## |
|---|
| 922 | #################### mrange ########################################## |
|---|
| 923 | # |
|---|
| 924 | # |
|---|
| 925 | # return a subset of selected columns, selected rows |
|---|
| 926 | # indexing always starts with 0. "end" can be used as an index. |
|---|
| 927 | # You are specifying reverse ordering for the selection when start>last. |
|---|
| 928 | # |
|---|
| 929 | proc mrange {m colstart collast {rowstart 0} {rowlast end}} { |
|---|
| 930 | set dim [dim_br m] |
|---|
| 931 | if { $dim < 2 } { promote_br m } |
|---|
| 932 | mrange_br m m $colstart $collast $rowstart $rowlast |
|---|
| 933 | return $m |
|---|
| 934 | } |
|---|
| 935 | |
|---|
| 936 | # |
|---|
| 937 | proc mrange_br {a_in c_out colstart collast {rowstart 0} {rowlast end}} { |
|---|
| 938 | upvar $a_in a |
|---|
| 939 | # matrix is assumed |
|---|
| 940 | set NR [lindex $a 1] |
|---|
| 941 | set NC [lindex $a 2] |
|---|
| 942 | foreach var {colstart collast} { |
|---|
| 943 | if { [set $var] == {end} } { set $var [expr {$NC - 1}] ; continue } |
|---|
| 944 | set $var [expr int([set $var])] |
|---|
| 945 | if { [set $var] < 0 } { error "column index should not be -ve"}\ |
|---|
| 946 | elseif { [set $var] >= $NC } { error "column index [set $var] > [expr {$NC - 1}]" } |
|---|
| 947 | } |
|---|
| 948 | if { $colstart <= $collast } { set colstep 1 } else {set colstep -1} |
|---|
| 949 | foreach var {rowstart rowlast} { |
|---|
| 950 | if { [set $var] == {end} } { set $var [expr {$NR - 1}] ; continue } |
|---|
| 951 | set $var [expr int([set $var])] |
|---|
| 952 | if { [set $var] < 0 } { error "row index should not be -ve"}\ |
|---|
| 953 | elseif { [set $var] >= $NR } { error "row index [set $var] > [expr {$NR - 1}]" } |
|---|
| 954 | } |
|---|
| 955 | if { $rowstart <= $rowlast } { set rowstep 1 } else {set rowstep -1} |
|---|
| 956 | set newNR [expr {$rowstep*($rowlast - $rowstart) + 1}] |
|---|
| 957 | set newNC [expr {$colstep*($collast - $colstart) + 1}] |
|---|
| 958 | set result [list 2 $newNR $newNC] |
|---|
| 959 | for {set row $rowstart} {1} {incr row $rowstep} { |
|---|
| 960 | for {set col $colstart} {1} {incr col $colstep} { |
|---|
| 961 | set index [expr {3 + $row*$NC + $col}] |
|---|
| 962 | lappend result [lindex $a $index] |
|---|
| 963 | if { $col == $collast } break |
|---|
| 964 | } |
|---|
| 965 | if { $row == $rowlast } break |
|---|
| 966 | } |
|---|
| 967 | upvar $c_out c |
|---|
| 968 | set c $result |
|---|
| 969 | return $c_out |
|---|
| 970 | } |
|---|
| 971 | |
|---|
| 972 | ############################ msolve ################################# |
|---|
| 973 | ############################ msolve ################################# |
|---|
| 974 | ############################ msolve ################################# |
|---|
| 975 | |
|---|
| 976 | # |
|---|
| 977 | # solve the matrix problem Ax = p for x, where p may be multiple columns |
|---|
| 978 | # |
|---|
| 979 | # when p is the identity matrix, the solution x, is the inverse of A |
|---|
| 980 | # |
|---|
| 981 | |
|---|
| 982 | |
|---|
| 983 | proc msolve {A p} { |
|---|
| 984 | promote_br p ;# a vector becomes Nx1 |
|---|
| 985 | |
|---|
| 986 | set N [lindex $A 2] |
|---|
| 987 | # paste the rhs columns on the end |
|---|
| 988 | join_cols_br A p |
|---|
| 989 | |
|---|
| 990 | msolve_br A |
|---|
| 991 | |
|---|
| 992 | # now peel off the solution which replaced the rhs columns |
|---|
| 993 | mrange_br A x $N end |
|---|
| 994 | return $x |
|---|
| 995 | } |
|---|
| 996 | |
|---|
| 997 | # |
|---|
| 998 | # Gauss elimination with partial pivoting |
|---|
| 999 | # |
|---|
| 1000 | proc msolve_br {A_in {tolerance 2.3e-16}} { |
|---|
| 1001 | upvar $A_in A |
|---|
| 1002 | set n [lindex $A 1] |
|---|
| 1003 | set NC [lindex $A 2] |
|---|
| 1004 | set p [expr {$NC - $n}] |
|---|
| 1005 | set Det 1.0 |
|---|
| 1006 | for {set j 0} {$j <= [expr {$n - 2}]} {incr j} { |
|---|
| 1007 | set j_plus1 [expr {$j + 1}] |
|---|
| 1008 | set jj [expr {3 + $j * $NC + $j}] |
|---|
| 1009 | set Ajj [lindex $A $jj] |
|---|
| 1010 | set s [expr {abs($Ajj)}] |
|---|
| 1011 | set k $j |
|---|
| 1012 | for {set h $j_plus1} {$h < $n} {incr h} { |
|---|
| 1013 | set Ahj [lindex $A [expr {3+ $h * $NC + $j}]] |
|---|
| 1014 | set Ahj [expr {abs($Ahj)}] |
|---|
| 1015 | if { $Ahj > $s } { |
|---|
| 1016 | set s $Ahj |
|---|
| 1017 | set k $h |
|---|
| 1018 | } |
|---|
| 1019 | } ;# end pivot search |
|---|
| 1020 | if { $k != $j } { ;# row interchange |
|---|
| 1021 | for {set i $j} {$i < $NC} {incr i} { |
|---|
| 1022 | set ki [expr {3+$k*$NC + $i}] |
|---|
| 1023 | set ji [expr {3+$j*$NC + $i}] |
|---|
| 1024 | set Aki [lindex $A $ki] |
|---|
| 1025 | set Aji [lindex $A $ji] |
|---|
| 1026 | lassign_br A $ki $Aji |
|---|
| 1027 | lassign_br A $ji $Aki |
|---|
| 1028 | } |
|---|
| 1029 | set Det [expr {0.0-$Det}] |
|---|
| 1030 | } |
|---|
| 1031 | set Ajj [lindex $A $jj] |
|---|
| 1032 | set Det [expr {$Det * $Ajj}] |
|---|
| 1033 | if { abs($Ajj) <= $tolerance } { |
|---|
| 1034 | error "matrix is computationally singular at a tolerance of $tolerance" |
|---|
| 1035 | } |
|---|
| 1036 | for {set k $j_plus1} {$k < $n} {incr k} { |
|---|
| 1037 | set kj [expr {3+$k*$NC + $j}] |
|---|
| 1038 | set Akj [lindex $A $kj] |
|---|
| 1039 | set Akj [expr {double($Akj)/$Ajj}] |
|---|
| 1040 | lassign_br A $kj $Akj ;# to form multiplier mkj |
|---|
| 1041 | for {set i $j_plus1} {$i < $NC} {incr i} { |
|---|
| 1042 | set ki [expr {3+$k*$NC + $i}] |
|---|
| 1043 | set Aki [lindex $A $ki] |
|---|
| 1044 | set ji [expr {3+$j*$NC + $i}] |
|---|
| 1045 | set Aji [lindex $A $ji] |
|---|
| 1046 | lassign_br A $ki [expr {$Aki - double($Akj)*$Aji}] |
|---|
| 1047 | } |
|---|
| 1048 | } ;# k |
|---|
| 1049 | } ;# j |
|---|
| 1050 | set nn [expr {3+($n-1)*$NC+($n-1)}] ;# n-1,n-1 |
|---|
| 1051 | set Ann [lindex $A $nn] |
|---|
| 1052 | set Det [expr {$Det * $Ann}] |
|---|
| 1053 | if { abs($Ann) < $tolerance } { |
|---|
| 1054 | error "matrix is computationally singular at a tolerance of $tolerance" |
|---|
| 1055 | } |
|---|
| 1056 | # factoring completed |
|---|
| 1057 | # mij stored in i>j |
|---|
| 1058 | # Rij stored in i<=j < n |
|---|
| 1059 | # fij stored in j= n ... n+p-1 |
|---|
| 1060 | # Back substitution for Rx=f |
|---|
| 1061 | for {set i $n} {$i < $NC} {incr i} { |
|---|
| 1062 | set ni [expr {3+($n-1)*$NC + $i}] ;# really n-1,i |
|---|
| 1063 | set Ani [lindex $A $ni] |
|---|
| 1064 | set Ani [expr {double($Ani)/$Ann}] |
|---|
| 1065 | lassign_br A $ni $Ani |
|---|
| 1066 | for {set j [expr {$n-2}]} {$j >= 0} {incr j -1} { |
|---|
| 1067 | set ji [expr {3+$j*$NC+$i}] |
|---|
| 1068 | set s [lindex $A $ji] |
|---|
| 1069 | for {set k [expr {$j+1}]} {$k < $n} {incr k} { |
|---|
| 1070 | set jk [expr {3+$j*$NC+$k}] |
|---|
| 1071 | set Ajk [lindex $A $jk] |
|---|
| 1072 | set ki [expr {3+$k*$NC+$i}] |
|---|
| 1073 | set Aki [lindex $A $ki] |
|---|
| 1074 | set s [expr {$s - $Ajk*$Aki}] |
|---|
| 1075 | } |
|---|
| 1076 | set jj [expr {3+$j*$NC+$j}] |
|---|
| 1077 | set Ajj [lindex $A $jj] |
|---|
| 1078 | lassign_br A $ji [expr {double($s)/$Ajj}] |
|---|
| 1079 | } |
|---|
| 1080 | } |
|---|
| 1081 | return $A_in |
|---|
| 1082 | } |
|---|
| 1083 | |
|---|
| 1084 | |
|---|
| 1085 | ########################## msum ##################################### |
|---|
| 1086 | ########################## msum ##################################### |
|---|
| 1087 | ########################## msum ##################################### |
|---|
| 1088 | # |
|---|
| 1089 | # compute the sums of each column, return a vector or scalar |
|---|
| 1090 | # call twice to get sum of columns and rows (set total [msum [msum $a]]) |
|---|
| 1091 | # |
|---|
| 1092 | # a is a vector or matrix |
|---|
| 1093 | # result is a vector or scalar |
|---|
| 1094 | proc msum a { |
|---|
| 1095 | if { [dim_br a] == 1 } { promote_br a } |
|---|
| 1096 | msum_br a sums |
|---|
| 1097 | # convert vector to scalar if only 1 column |
|---|
| 1098 | return [demote $sums] |
|---|
| 1099 | } |
|---|
| 1100 | |
|---|
| 1101 | # a_in is matrix |
|---|
| 1102 | # sums_out writes a vector |
|---|
| 1103 | proc msum_br {a_in sums_out} { |
|---|
| 1104 | upvar $a_in a |
|---|
| 1105 | if { [dim_br a] != 2 } { error "expecting matrix" } |
|---|
| 1106 | set NR [lindex $a 1] |
|---|
| 1107 | set NC [lindex $a 2] |
|---|
| 1108 | # vector results are returned |
|---|
| 1109 | set sums [list 2 $NC 0] |
|---|
| 1110 | for {set j 0} {$j < $NC} {incr j} { |
|---|
| 1111 | set sum 0.0 |
|---|
| 1112 | set index [expr {3 + $j}] |
|---|
| 1113 | for {set i 0} {$i < $NR} {incr i} { |
|---|
| 1114 | set Aij [lindex $a $index] |
|---|
| 1115 | set sum [expr {$sum + $Aij}] |
|---|
| 1116 | incr index $NC |
|---|
| 1117 | } |
|---|
| 1118 | lappend sums $sum |
|---|
| 1119 | } |
|---|
| 1120 | upvar $sums_out m |
|---|
| 1121 | set m $sums |
|---|
| 1122 | return $sums_out |
|---|
| 1123 | } |
|---|
| 1124 | |
|---|
| 1125 | |
|---|
| 1126 | ################################ msvd ########################## |
|---|
| 1127 | ################################ msvd ########################## |
|---|
| 1128 | ################################ msvd ########################## |
|---|
| 1129 | # |
|---|
| 1130 | # Singular Value Decomposition |
|---|
| 1131 | # decompose matrix A into (U)(S)(Vtrans) where |
|---|
| 1132 | # A[m,n] is the original matrix |
|---|
| 1133 | # U[m,n] has orthogonal columns (Ut)(U) = (1(k) |
|---|
| 1134 | # and multiplies to an identity matrix ... |
|---|
| 1135 | # supplemented with zeroes if needed 0(n-k)) |
|---|
| 1136 | # V[n,n] is orthogonal (V)(Vtran) = [mident $n] |
|---|
| 1137 | # the eigenvectors representing the principal components |
|---|
| 1138 | # S is diagonal with the positive singular values of A |
|---|
| 1139 | # |
|---|
| 1140 | # when you have V,S,U |
|---|
| 1141 | # A[m,n]V[n,n] = B[m,n] is a transformation of A to orthogonal columns, B |
|---|
| 1142 | # B[m,n] = U[m,n]S[n,n] |
|---|
| 1143 | # square S and divide by (m-1) to get PCA eigenvalues |
|---|
| 1144 | # |
|---|
| 1145 | proc msvd A { |
|---|
| 1146 | msvd_br A S V |
|---|
| 1147 | puts U:\n[show $A %12.4f] |
|---|
| 1148 | puts \nS:\n[show $S %12.4f] |
|---|
| 1149 | puts \nV:\n[show $V %12.4f] |
|---|
| 1150 | return $V |
|---|
| 1151 | } |
|---|
| 1152 | # |
|---|
| 1153 | proc msvd_br {A_in_U_out S_out V_out {epsilon 2.3e-16}} { |
|---|
| 1154 | upvar $A_in_U_out A |
|---|
| 1155 | set m [lindex $A 1] |
|---|
| 1156 | set n [lindex $A 2] |
|---|
| 1157 | set tolerance [expr {$epsilon * $epsilon* $m * $n}] |
|---|
| 1158 | upvar $V_out V |
|---|
| 1159 | set V [mident $n] |
|---|
| 1160 | upvar $S_out z |
|---|
| 1161 | |
|---|
| 1162 | # top of the iteration |
|---|
| 1163 | set count 1 |
|---|
| 1164 | for {set isweep 0} {$isweep < 30 && $count > 0} {incr isweep} { |
|---|
| 1165 | set count [expr {$n*($n-1)/2}] ;# count of rotations in a sweep |
|---|
| 1166 | for {set j 0} {$j < [expr {$n-1}]} {incr j} { |
|---|
| 1167 | for {set k [expr {$j+1}]} {$k < $n} {incr k} { |
|---|
| 1168 | set p [set q [set r 0.0]] |
|---|
| 1169 | for {set i 0} {$i < $m} {incr i} { |
|---|
| 1170 | set ij [expr {3+$i*$n+$j}] |
|---|
| 1171 | set ik [expr {3+$i*$n+$k}] |
|---|
| 1172 | set Aij [lindex $A $ij] |
|---|
| 1173 | set Aik [lindex $A $ik] |
|---|
| 1174 | set p [expr {$p + $Aij*$Aik}] |
|---|
| 1175 | set q [expr {$q + $Aij*$Aij}] |
|---|
| 1176 | set r [expr {$r + $Aik*$Aik}] |
|---|
| 1177 | } |
|---|
| 1178 | if { $q < $r } { |
|---|
| 1179 | set c 0.0 |
|---|
| 1180 | set s 1.0 |
|---|
| 1181 | }\ |
|---|
| 1182 | elseif { $q * $r == 0.0 } { ;# underflow of small elements |
|---|
| 1183 | incr count -1 |
|---|
| 1184 | continue |
|---|
| 1185 | }\ |
|---|
| 1186 | elseif { ($p*$p)/($q*$r) < $tolerance } { ;# cols j,k are orthogonal |
|---|
| 1187 | incr count -1 |
|---|
| 1188 | continue |
|---|
| 1189 | }\ |
|---|
| 1190 | else { |
|---|
| 1191 | set q [expr {$q-$r}] |
|---|
| 1192 | set v [expr {sqrt(4.0*$p*$p + $q*$q)}] |
|---|
| 1193 | set c [expr {sqrt(($v+$q)/(2.0*$v))}] |
|---|
| 1194 | set s [expr {$p/($v*$c)}] |
|---|
| 1195 | # s == sine of rotation angle, c == cosine |
|---|
| 1196 | } |
|---|
| 1197 | # rotation of A |
|---|
| 1198 | for {set i 0} {$i < $m} {incr i} { |
|---|
| 1199 | set ij [expr {3+$i*$n+$j}] |
|---|
| 1200 | set ik [expr {3+$i*$n+$k}] |
|---|
| 1201 | set Aij [lindex $A $ij] |
|---|
| 1202 | set Aik [lindex $A $ik] |
|---|
| 1203 | lassign_br A $ij [expr {$Aij*$c + $Aik*$s}] |
|---|
| 1204 | lassign_br A $ik [expr {$Aik*$c - $Aij*$s}] |
|---|
| 1205 | } |
|---|
| 1206 | # rotation of V |
|---|
| 1207 | for {set i 0} {$i < $n} {incr i} { |
|---|
| 1208 | set ij [expr {3+$i*$n+$j}] |
|---|
| 1209 | set ik [expr {3+$i*$n+$k}] |
|---|
| 1210 | set Vij [lindex $V $ij] |
|---|
| 1211 | set Vik [lindex $V $ik] |
|---|
| 1212 | lassign_br V $ij [expr {$Vij*$c + $Vik*$s}] |
|---|
| 1213 | lassign_br V $ik [expr {$Vik*$c - $Vij*$s}] |
|---|
| 1214 | } |
|---|
| 1215 | } ;# k |
|---|
| 1216 | } ;# j |
|---|
| 1217 | #puts "pass=$isweep skipped rotations=$count" |
|---|
| 1218 | } ;# isweep |
|---|
| 1219 | |
|---|
| 1220 | set z [list 2 $n 0] |
|---|
| 1221 | for {set j 0} {$j < $n} {incr j} { |
|---|
| 1222 | set q 0.0 |
|---|
| 1223 | for {set i 0} {$i < $m} {incr i} { |
|---|
| 1224 | set ij [expr {3+$i*$n+$j}] |
|---|
| 1225 | set Aij [lindex $A $ij] |
|---|
| 1226 | set q [expr {$q+$Aij*$Aij}] |
|---|
| 1227 | } |
|---|
| 1228 | set q [expr {sqrt($q)}] |
|---|
| 1229 | lappend z $q |
|---|
| 1230 | if { $q >= $tolerance } { |
|---|
| 1231 | for {set i 0} {$i < $m} {incr i} { |
|---|
| 1232 | set ij [expr {3+$i*$n+$j}] |
|---|
| 1233 | set Aij [lindex $A $ij] |
|---|
| 1234 | lassign_br A $ij [expr {$Aij/$q}] |
|---|
| 1235 | } |
|---|
| 1236 | } |
|---|
| 1237 | } ;# j |
|---|
| 1238 | return [list $A_in_U_out $S_out $V_out] |
|---|
| 1239 | } |
|---|
| 1240 | |
|---|
| 1241 | |
|---|
| 1242 | ######### mat_unary_op, madjust, moffset, mscale #################### |
|---|
| 1243 | ######### mat_unary_op, madjust, moffset, mscale #################### |
|---|
| 1244 | ######### mat_unary_op, madjust, moffset, mscale #################### |
|---|
| 1245 | # |
|---|
| 1246 | # matrix unary operations like scaling |
|---|
| 1247 | # |
|---|
| 1248 | proc mat_unary_op {a op} { |
|---|
| 1249 | set dimA [dim_br a] |
|---|
| 1250 | if { $dimA > 2} { |
|---|
| 1251 | error "cannot handle > 2D data" |
|---|
| 1252 | } |
|---|
| 1253 | if { $dimA < 1 } { promote_br a } |
|---|
| 1254 | mat_unary_op_br a $op c |
|---|
| 1255 | return $c |
|---|
| 1256 | } |
|---|
| 1257 | |
|---|
| 1258 | proc moffset {a delta} {return [mat_unary_op $a "expr $delta +"] } |
|---|
| 1259 | proc mscale {a delta} { return [mat_unary_op $a "expr $delta *"] } |
|---|
| 1260 | proc madjust {a scale offset} { return [mat_unary_op $a "expr $offset + $scale *"] } |
|---|
| 1261 | proc mround_ij {eps aij} { ;# round off a number if close to an integer |
|---|
| 1262 | if { abs($aij - int($aij)) <= $eps } { return [expr {double(int($aij))}] } |
|---|
| 1263 | return $aij |
|---|
| 1264 | } |
|---|
| 1265 | proc mround {a {eps 1.0e-8}} { return [mat_unary_op $a [list mround_ij $eps]] } |
|---|
| 1266 | |
|---|
| 1267 | # perform a unary operation on the elements of a vector or matrix |
|---|
| 1268 | # the value is lappended to your op and the result is executed |
|---|
| 1269 | # "expr 0.5 *" divides by 2 |
|---|
| 1270 | # you can use a procedure for more complex formulae |
|---|
| 1271 | # |
|---|
| 1272 | proc mat_unary_op_br {a_in op {c_out {}}} { |
|---|
| 1273 | upvar $a_in a |
|---|
| 1274 | if { $c_out == {}} { set c_out $a_in } |
|---|
| 1275 | upvar $c_out c |
|---|
| 1276 | set NR [lindex $a 1] |
|---|
| 1277 | set MC [lindex $a 2] |
|---|
| 1278 | set result [list 2 $NR $MC] |
|---|
| 1279 | set imax [llength $a] |
|---|
| 1280 | for {set i 3} {$i < $imax} {incr i} { |
|---|
| 1281 | lappend result [eval [concat $op [lindex $a $i]]] |
|---|
| 1282 | } |
|---|
| 1283 | set c $result ;# only changes a_in,c if no error |
|---|
| 1284 | return $c_out |
|---|
| 1285 | } |
|---|
| 1286 | ########################### promote ################################# |
|---|
| 1287 | ########################### promote ################################# |
|---|
| 1288 | ########################### promote ################################# |
|---|
| 1289 | # |
|---|
| 1290 | # promote a scalar or vector to an array |
|---|
| 1291 | # vec[N] promoted to Nx1 array |
|---|
| 1292 | # |
|---|
| 1293 | proc promote x { |
|---|
| 1294 | promote_br x |
|---|
| 1295 | return $x |
|---|
| 1296 | } |
|---|
| 1297 | |
|---|
| 1298 | proc promote_br {name_in {name_out {}}} { |
|---|
| 1299 | upvar $name_in x |
|---|
| 1300 | if { $name_out == {}} {set name_out $name_in } |
|---|
| 1301 | upvar $name_out y |
|---|
| 1302 | set dim [dim_br x] |
|---|
| 1303 | if { $dim == 0 } { |
|---|
| 1304 | set y [list 2 1 1 $x] |
|---|
| 1305 | return $name_out |
|---|
| 1306 | } |
|---|
| 1307 | if { $dim == 1 } { ;# 2 n 0 v1 v2 ... |
|---|
| 1308 | set n [lindex $x 1] |
|---|
| 1309 | if { $name_out != $name_in } { set y [lreplace $x 2 2 1] }\ |
|---|
| 1310 | else { lassign_br y 2 1 } |
|---|
| 1311 | return $name_out |
|---|
| 1312 | } |
|---|
| 1313 | if { $name_out != $name_in } {set y $x} |
|---|
| 1314 | return $name_out |
|---|
| 1315 | } |
|---|
| 1316 | |
|---|
| 1317 | |
|---|
| 1318 | # demote an Nx1 or 1XN matrix to a vect[n] |
|---|
| 1319 | # demote a vect[1] to a scalar |
|---|
| 1320 | # call twice to demote 1x1 matrix to scalar |
|---|
| 1321 | proc demote x { |
|---|
| 1322 | demote_br x |
|---|
| 1323 | return $x |
|---|
| 1324 | } |
|---|
| 1325 | |
|---|
| 1326 | proc demote_br {a_in {c_out {}}} { |
|---|
| 1327 | upvar $a_in a |
|---|
| 1328 | if { $c_out == {}} { set c_out $a_in } |
|---|
| 1329 | upvar $c_out c |
|---|
| 1330 | if { $c_out != $a_in } { set c $a } |
|---|
| 1331 | # demote 1XN or Nx1 to vector |
|---|
| 1332 | set dim [dim_br a] |
|---|
| 1333 | if { $dim == 0 || $dim > 2 } {return $c_out } |
|---|
| 1334 | set nr [lindex $a 1] |
|---|
| 1335 | if { $dim == 2 } { |
|---|
| 1336 | set mc [lindex $a 2] |
|---|
| 1337 | if { $mc == 1 } { |
|---|
| 1338 | set c [lreplace $a 2 2 0] |
|---|
| 1339 | }\ |
|---|
| 1340 | elseif { $nr == 1 } { |
|---|
| 1341 | set c [lreplace $a 1 2 $mc 0] |
|---|
| 1342 | } |
|---|
| 1343 | return $c_out |
|---|
| 1344 | } |
|---|
| 1345 | # demote vector1 to scalar |
|---|
| 1346 | if { $dim == 1 } { |
|---|
| 1347 | if { $nr == 1 } { |
|---|
| 1348 | set c [lindex $a 3] |
|---|
| 1349 | } |
|---|
| 1350 | } |
|---|
| 1351 | return $c_out |
|---|
| 1352 | } |
|---|
| 1353 | |
|---|
| 1354 | ############################ show ################################### |
|---|
| 1355 | ############################ show ################################### |
|---|
| 1356 | ############################ show ################################### |
|---|
| 1357 | |
|---|
| 1358 | # |
|---|
| 1359 | # Return a formatted string representation for an lalf operand |
|---|
| 1360 | # the format argument is optionally used per-element |
|---|
| 1361 | # with the Tcl format command. Examples: |
|---|
| 1362 | # %.4f gives 4 decimal digit fixed point |
|---|
| 1363 | # %12.4e provides 4 decimal digit exponential format with a fixed |
|---|
| 1364 | # width |
|---|
| 1365 | # |
|---|
| 1366 | # The col_join argument can be used to separate values with commas |
|---|
| 1367 | # and tabs, etc. Similarly the row_join allows you to define the |
|---|
| 1368 | # row separation string. |
|---|
| 1369 | # |
|---|
| 1370 | proc show {x {format {}} {col_join { }} {row_join \n}} { |
|---|
| 1371 | show_br x x $format $col_join $row_join |
|---|
| 1372 | return $x |
|---|
| 1373 | } |
|---|
| 1374 | |
|---|
| 1375 | proc show_br {name_in {name_out {}} {format {}} {col_join { }} {row_join \n}} { |
|---|
| 1376 | upvar $name_in x |
|---|
| 1377 | if { $name_out == {}} { set name_out $name_in } |
|---|
| 1378 | upvar $name_out result |
|---|
| 1379 | set dim [dim_br x] |
|---|
| 1380 | if { $dim == 0 } { |
|---|
| 1381 | if { $format != {} } {set result [format $format $x] }\ |
|---|
| 1382 | else { set result $x } |
|---|
| 1383 | return $name_out} |
|---|
| 1384 | if { $dim == 1 } { ;# 2 n 0 v1 v2 ... |
|---|
| 1385 | if { $format == {}} { |
|---|
| 1386 | set result [join [lrange $x 3 end] $col_join] |
|---|
| 1387 | }\ |
|---|
| 1388 | else { |
|---|
| 1389 | set temp [list] |
|---|
| 1390 | foreach item [lrange $x 3 end] { |
|---|
| 1391 | lappend temp [format $format $item] |
|---|
| 1392 | } |
|---|
| 1393 | set result [join $temp $col_join] |
|---|
| 1394 | } |
|---|
| 1395 | return $name_out |
|---|
| 1396 | } |
|---|
| 1397 | if { $dim == 2} { |
|---|
| 1398 | if { $format != {}} { |
|---|
| 1399 | mat_unary_op_br x [list format $format] result |
|---|
| 1400 | } \ |
|---|
| 1401 | else { |
|---|
| 1402 | set result $x |
|---|
| 1403 | } |
|---|
| 1404 | set rows [list] |
|---|
| 1405 | set NR [lindex $result 1] |
|---|
| 1406 | set MC [lindex $result 2] |
|---|
| 1407 | for {set row 0} {$row < $NR} {incr row} { |
|---|
| 1408 | set start [expr {3 + $row*$MC}] |
|---|
| 1409 | set end [expr {$start + $MC - 1}] |
|---|
| 1410 | lappend rows [join [lrange $result $start $end] $col_join] |
|---|
| 1411 | } |
|---|
| 1412 | set result [join $rows $row_join] |
|---|
| 1413 | return $name_out |
|---|
| 1414 | } |
|---|
| 1415 | |
|---|
| 1416 | error "> 2-D TBD" |
|---|
| 1417 | } |
|---|
| 1418 | |
|---|
| 1419 | |
|---|
| 1420 | ########################### transpose ############################### |
|---|
| 1421 | ########################### transpose ############################### |
|---|
| 1422 | ########################### transpose ############################### |
|---|
| 1423 | # |
|---|
| 1424 | # exchange [i,j] with [j,i] |
|---|
| 1425 | # |
|---|
| 1426 | # a vector is promoted to a 1xN array by transpose (1 row x N col) |
|---|
| 1427 | # |
|---|
| 1428 | #% set a {2 2 3 01 02 03 11 12 13} |
|---|
| 1429 | #% show $a |
|---|
| 1430 | # 01 02 03 |
|---|
| 1431 | # 11 12 13 |
|---|
| 1432 | #% transpose $a |
|---|
| 1433 | # 2 3 2 01 11 02 12 03 13 |
|---|
| 1434 | #% show [transpose $a] |
|---|
| 1435 | # 01 11 |
|---|
| 1436 | # 02 12 |
|---|
| 1437 | # 03 13 |
|---|
| 1438 | # |
|---|
| 1439 | proc transpose x { ;# call by value |
|---|
| 1440 | transpose_br x |
|---|
| 1441 | return $x |
|---|
| 1442 | } |
|---|
| 1443 | |
|---|
| 1444 | # transpose |
|---|
| 1445 | # call by reference - default output is to overwrite input data |
|---|
| 1446 | # |
|---|
| 1447 | proc transpose_br {name_in {name_out {}}} { |
|---|
| 1448 | upvar $name_in x |
|---|
| 1449 | if { $name_out == {}} {set name_out $name_in } |
|---|
| 1450 | upvar $name_out y |
|---|
| 1451 | set dim [dim_br x] |
|---|
| 1452 | if { $dim == 0 } { return } |
|---|
| 1453 | if { $dim == 1 } { ;# {2 N 0 ...} vector to 1xN |
|---|
| 1454 | set N [lindex $x 1] |
|---|
| 1455 | set y [lreplace $x 1 2 1 $N] |
|---|
| 1456 | return $name_out |
|---|
| 1457 | } |
|---|
| 1458 | if { $dim == 2 } { |
|---|
| 1459 | |
|---|
| 1460 | set NR [lindex $x 1] |
|---|
| 1461 | set MC [lindex $x 2] |
|---|
| 1462 | set result [list 2 $MC $NR] |
|---|
| 1463 | for {set j 0} {$j < $MC} {incr j} { ;# j == col in source |
|---|
| 1464 | for { set i 0 } { $i < $NR } { incr i } { ;# i == row in source |
|---|
| 1465 | lappend result [lindex $x [expr {3 + $i * $MC + $j}]] |
|---|
| 1466 | } |
|---|
| 1467 | } |
|---|
| 1468 | set y $result |
|---|
| 1469 | return $name_out |
|---|
| 1470 | } |
|---|
| 1471 | error "transpose cannot handle >2D data" |
|---|
| 1472 | } |
|---|
| 1473 | |
|---|
| 1474 | ########################### vdiag ############################# |
|---|
| 1475 | ########################### vdiag ############################# |
|---|
| 1476 | ########################### vdiag ############################# |
|---|
| 1477 | # |
|---|
| 1478 | # Create a vector from the diagonal elements of a matrix |
|---|
| 1479 | # |
|---|
| 1480 | proc vdiag m { |
|---|
| 1481 | promote_br m |
|---|
| 1482 | vdiag_br m v |
|---|
| 1483 | return $v |
|---|
| 1484 | } |
|---|
| 1485 | |
|---|
| 1486 | proc vdiag_br {name_m_in name_v_out} { |
|---|
| 1487 | upvar $name_m_in m |
|---|
| 1488 | upvar $name_v_out v |
|---|
| 1489 | set N [lindex $m 1] |
|---|
| 1490 | set C [lindex $m 2] |
|---|
| 1491 | if { $C < $N } { set N $C } |
|---|
| 1492 | set v [list 2 $N 0] |
|---|
| 1493 | for {set i 0} {$i < $N} {incr i} { |
|---|
| 1494 | set ij [expr {$i*($C + 1) + 3}] |
|---|
| 1495 | lappend v [lindex $m $ij] |
|---|
| 1496 | } |
|---|
| 1497 | |
|---|
| 1498 | return $name_v_out |
|---|
| 1499 | } |
|---|
| 1500 | |
|---|
| 1501 | |
|---|
| 1502 | ########################### vtrim ############################# |
|---|
| 1503 | ########################### vtrim ############################# |
|---|
| 1504 | ########################### vtrim ############################# |
|---|
| 1505 | |
|---|
| 1506 | # |
|---|
| 1507 | # for a vector, just return the actual vector elements |
|---|
| 1508 | # trim away the dimension and size data in the front |
|---|
| 1509 | # |
|---|
| 1510 | proc vtrim x { |
|---|
| 1511 | vtrim_br x |
|---|
| 1512 | return $x |
|---|
| 1513 | } |
|---|
| 1514 | |
|---|
| 1515 | proc vtrim_br {name_in {name_out {}}} { |
|---|
| 1516 | upvar $name_in x |
|---|
| 1517 | if { $name_out == {}} {set name_out $name_in } |
|---|
| 1518 | upvar $name_out y |
|---|
| 1519 | set y [lrange $x 3 end] |
|---|
| 1520 | return $name_out |
|---|
| 1521 | } |
|---|
| 1522 | |
|---|
| 1523 | |
|---|
| 1524 | |
|---|
| 1525 | |
|---|
| 1526 | |
|---|
| 1527 | |
|---|
| 1528 | |
|---|
| 1529 | |
|---|
| 1530 | |
|---|
| 1531 | |
|---|
| 1532 | |
|---|
| 1533 | } ;# end of namespace |
|---|
