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 |
---|