source: trunk/la1.0/la.tcl @ 931

Last change on this file since 931 was 931, checked in by toby, 11 years ago

include rest of files

File size: 50.0 KB
Line 
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
18package 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.
153global La_namespace
154if { ![info exists La_namespace] } { set La_namespace La }
155
156namespace eval $La_namespace {
157
158  namespace export \
159dim dim_br \
160dotprod dotprod_br \
161join_cols join_cols_br \
162join_rows join_rows_br\
163lassign lassign_br \
164madd mdiv mprod msub mat_binary_op mat_binary_op_br \
165mathprec \
166mcols mcols_br\
167mdiag \
168mdingdong \
169mevsvd mevsvd_br \
170mhilbert \
171mident \
172mlssvd \
173mmult mmult_br \
174mnorms mnorms_br mnormalize mnormalize_br \
175mrange mrange_br \
176mround mround_ij \
177mrows mrows_br\
178msolve msolve_br \
179msum msum_br \
180msvd msvd_br\
181mat_unary_op madjust moffset mscale mat_unary_op_br \
182promote promote_br demote demote_br\
183show show_br \
184transpose transpose_br \
185vdiag vdiag_br \
186vtrim 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
207proc dim x {
208    return [dim_br x]
209    }
210
211proc 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#
249proc 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#
256proc 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#
285proc 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   
298proc 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#
332proc 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   
345proc 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#
371proc 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)
377proc 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.
384if { [info commands ${La_namespace}::lassign_br] == {}} {
385proc 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#
399proc 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
412proc madd {a b} { return [mat_binary_op $a $b +] }
413proc mdiv {a b} { return [mat_binary_op $a $b {*1.0/}] }
414proc mprod {a b} { return [mat_binary_op $a $b *] }
415proc 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#
425proc 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#
463proc 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#
487proc mcols m { return [lindex $m 2] }
488proc mcols_br name_m { upvar $name_m m ; return [lindex $m 2] }
489proc mrows m { return [lindex $m 1] } 
490proc 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#
498proc 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#
522proc 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#
540proc mevsvd {A {eps 2.3e-16}} {
541    mevsvd_br A evals
542puts "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
550proc 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#
664proc 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#
681proc 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#
709proc 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#   
770proc 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
783proc 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#
818proc mnorms a {
819    mnorms_br a means sigmas
820    return [transpose [join_cols $means $sigmas]]
821    }
822
823proc 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   
850proc 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#
859proc 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
889if { [info commands mean] != "mean" } {
890proc ::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}
902if { [info commands stddev] != "stddev" } {
903proc ::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#
929proc 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#
937proc 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
983proc 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#
1000proc 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
1094proc 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
1103proc 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#
1145proc msvd A {
1146    msvd_br A S V
1147puts U:\n[show $A %12.4f] 
1148puts \nS:\n[show $S %12.4f]
1149puts \nV:\n[show $V %12.4f]
1150    return $V
1151    }
1152#
1153proc 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#
1248proc 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
1258proc moffset {a delta} {return [mat_unary_op $a "expr $delta +"] }
1259proc mscale {a delta} { return [mat_unary_op $a "expr $delta *"] }
1260proc madjust {a scale offset} { return [mat_unary_op $a "expr $offset + $scale *"] }
1261proc 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    }
1265proc 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#
1272proc 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#
1293proc promote x {
1294    promote_br x
1295    return $x
1296    }
1297
1298proc 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
1321proc demote x {
1322    demote_br x
1323    return $x
1324    }
1325
1326proc 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#
1370proc show {x {format {}} {col_join { }} {row_join \n}} {
1371    show_br x x $format $col_join $row_join
1372    return $x
1373    }
1374
1375proc 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#
1439proc 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#
1447proc 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#
1480proc vdiag m {
1481    promote_br m
1482    vdiag_br m v
1483    return $v
1484    }
1485
1486proc 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#
1510proc vtrim x {
1511    vtrim_br x
1512    return $x
1513    }
1514
1515proc 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
Note: See TracBrowser for help on using the repository browser.