Version 9 of Numerical array operations

Updated 2004-02-27 17:09:15

if { 0 } { Arjen Markus (18 february 2004) Several languages, among which Fortran 90, have array operations, to use Fortran 90 syntax:

   real, dimension(1:100) :: a, b, c

   b = 1.0
   do i = 1,size(a)
      a(i) = i
   enddo

   c= b / a 

would fill the array b with the value 1.0, in the do-loop each element of a is given a value (otherwise the last assignment is pretty dull) and each element of the array c is assigned the reciprocal value of the corresponding element of a.

In other words: arrays need not be manipulated only through do-loops (or for-loops, if you like).

This intrigued me: can we get the same thing into Tcl? The answer is obvious: yes!

The script below does very little error handling and it is currently limited to one-dimensional arrays (i.e., Tcl lists of simple numbers) only. But it is just a proof of concept :)


Note that arrays here are not Tcl's hash tables, but numerically indexed collections - for which one best uses lists in Tcl.

AM It is not always easy to come up with words that have an unambiguous meaning. I did not want to use "list" as that would have reminded people too much of other list operations...

See also the NAP package and the LA package (to a lesser extent).

AM (24 february 2004) The next steps in this experiment:

  • Add more operations (I have already done a few)
  • Measuring the performance (I have already done a few measurements)
  • Make it possible to use multi-dimensional arrays - I intend to use both the approach taken by LA and the one in Playing APL

}

  # numarray.tcl --
 #    Array operations on lists of numbers
 #
 namespace eval ::numarray {
    variable __v__
 }

 proc ::numarray::numrange { varname vmin vmax vstep } {
    upvar $varname var

    if { $vmax < $vmin } {
       return -code error "Minimum must be smaller than maximum"
    }
    if { $vstep <= 0.0 } {
       return -code error "Step must be positive"
    }

    set value $vmin
    set var   {}
    while { $value < $vmax+0.5*$vstep } {
       lappend var $value
       set value [expr {$value+$vstep}]
    }
 }

 proc ::numarray::numset { varname expression } {
    upvar $varname var

    #
    # Replace the substrings "$var" by a corresponding
    # list value - if the variable is a list
    #
    set _length_ 1

    while { [string first \$ $expression] >= 0 } {
       regexp {\$(\w+)} $expression ==> _v_
       catch { upvar $_v_ $_v_ }
       if { [llength [set $_v_]] > 1 } {
          set _length_ [llength [set $_v_]]
          regsub {\$\w+} $expression "\[lindex \@$_v_ @_i_\]" expression
       } else {
          regsub {\$\w+} $expression "@$_v_" expression
       }
    }

    set expression [string map {@ $} $expression]

    set _result_ {}
    for { set _i_ 0 } { $_i_ < $_length_ } { incr _i_ } {
       lappend _result_ [expr $expression]
    }
    set var $_result_
 }

 #
 # Using numset to create other array operations:
 # efficient coding, but there is performance to be gained
 # by reimplementing them
 #
 proc ::numarray::all {condition} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $condition]

    puts "All: $__v__"
    expr {[lsearch $__v__ 0] < 0}
 }

 proc ::numarray::any {condition} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $condition]

    expr {[lsearch $__v__ 1] >= 0}
 }

 proc ::numarray::count {condition} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $condition]

    regexp -all {1} $__v__
 }

 proc ::numarray::maxval {expression} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $expression]

    set max [lindex $__v__ 0]
    foreach v $__v__ {
       if { $max < $v } {
          set max $v
       }
    }
    return $max
 }

 proc ::numarray::minval {expression} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $expression]

    set min [lindex $__v__ 0]
    foreach v $__v__ {
       if { $min > $v } {
          set min $v
       }
    }
    return $min
 }

 proc ::numarray::maxloc {expression} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $expression]

    set max    [lindex $__v__ 0]
    set maxidx 0
    set idx    0
    foreach v $__v__ {
       if { $max < $v } {
          set max    $v
          set maxidx $idx
       }
       incr idx
    }
    return $maxidx
 }

 proc ::numarray::minloc {expression} {
    variable __v__
    uplevel [list ::numarray::numset ::numarray::__v__ $expression]

    set min    [lindex $__v__ 0]
    set minidx 0
    set idx    0
    foreach v $__v__ {
       if { $min > $v } {
          set min    $v
          set minidx $idx
       }
       incr idx
    }
    return $minidx
 }

 #
 # A few tests
 #
 catch { console show }

 ::numarray::numrange v 0.0 10.0 1.1
 puts "v: $v"
 ::numarray::numset w {$v*$v}

 set p 1.0
 ::numarray::numset v {$v+$w+$p}

 puts "v: $v"
 puts "w: $w"

 ::numarray::numrange v 0.0 10.0 2.5
 ::numarray::numrange w 1.0 10.0 2.0
 puts "Examine the values v = $v"
 puts "All values > -1? [::numarray::all {$v > -1}]"
 puts "All values < 4? [::numarray::all {$v < 4}]"
 puts "Any values < 4? [::numarray::any {$v < 4}]"
 puts "How many values < 4? [::numarray::count {$v < 4}]"
 puts "Now also: w = $w"
 puts "Any values of v < the corresponding values of w? [::numarray::any {$v < $w}]"

 puts "Maximum and minimum of v*w: [::numarray::maxval {$v*$w}], [::numarray::minval {$v*$w}]"
 puts "Maximum and minimum of v*w - where: [::numarray::maxloc {$v*$w}], [::numarray::minloc {$v*$w}]"

 #
 # Simple measurements
 # Note: the foreach loop must be in a procedure, otherwise it does not
 # get compiled and the measurement is biased. (Thanks to [MS])
 #
 proc testloop { v w } { 
    set r {}
    foreach a $v b $w {
       lappend r [expr {$a+$b+$a*$b}]
    }
    return $r
 }

 ::numarray::numrange v 0.0 10.0 0.001
 ::numarray::numrange w 1.0 11.0 0.001
 puts [time {::numarray::numset a {$v+$w+$v*$w}} 10]
 puts [time {testloop $v $w} 100]

 ::numarray::numrange v 0.0 10.0 0.01
 ::numarray::numrange w 1.0 11.0 0.01
 puts [time {::numarray::numset a {$v+$w+$v*$w}} 100]
 puts [time {testloop $v $w} 100]

 ::numarray::numrange v 0.0 10.0 0.1
 ::numarray::numrange w 1.0 11.0 0.1
 puts [time {::numarray::numset a {$v+$w+$v*$w}} 1000]
 puts [time {testloop $v $w} 1000]

 #
 # Displaying a curve ...
 #
 catch {

 proc merge { list1 list2 } {
    set result {}
    foreach e1 $list1 e2 $list2 {
       lappend result $e1 $e2
    }
    return $result
 }

 package require Tk

 ::numarray::numrange phi 0.0 [expr {2.0*3.1415926}] [expr {0.02*3.1415926}]
 ::numarray::numset   rad {1.0+cos($phi)}
 ::numarray::numset   x   {100 + 100 * $rad * cos($phi)}
 ::numarray::numset   y   {200 + 100 * $rad * sin($phi)}

 pack [canvas .c -bg white -width 400 -height 400]
 .c create line [merge $x $y]

 }

RS: Come to think, multi-dimensional "arrays" are best represented with nested lists, the same as trees, so in Tcl, the concepts of "matrix" and "tree" converge... and could both be handled by either recursive iteration, or a flat foreach, in conjunction with lset and lindex, over the index vectors...

AM IIRC, Ed Hume chose for a sort of tagged lists for performance reasons. Also, the method he chose makes it easy and natural to represent column and row vectors.

The method:

   {2 3 2 1.0 0.0 0.0 2.0 3.0 4.0}

represents a two-dimensional array (the first element is the number of dimensions) of 3 columns and 2 rows (the next elements give the array dimensions). So:

  / 1.0  0.0  0.0 \
  \ 2.0  3.0  4.0 /

The method can be extended to any number of dimensions you like and such a structure makes element-wise operations easy and fast, whereas accessing an individual element is also very easy.


See also lexpr


[ Category Concept

Category Language

]