Numerical array operations

Arjen Markus (18 february 2004) Several languages, among which Fortran 90/95, 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

TV (Mar 5, 2004) What's the actual problem?

 /* matrix.c */
 #include<stdio.h>

 #define MAX 10000000
 #define nopr

 double a[MAX], b[MAX], c[MAX];
 int i,j, l=MAX;
 #define FILLV(m,l,v) { for (i=0; i<l; i++) m[i] = v; };
 #define DIVV(m,l,n,o) { for (i=0; i<l; i++) m[i] = n[i]/o[i]; };

 main()
 {
   FILLV(b,l,1);
   FILLV(a,l,(double) i+1);
   DIVV(c,l,b,a);
 #ifdef pr
   for (j=0; j<l; j++)
      printf("%d %lf\n", j, c[j]);
 #endif  /* pr */
 }

Under unix/linux(/OSX?)/cygwin (on windows XX) compile with:

   gcc -o matrix matrix.c -lm

Run:

   time matrix

Result:

 real    0m0.900s
 user    0m0.370s
 sys     0m0.530s

900,000/10,000,000 = 0.09 USec/iter

reducing to floats makes for another halving of time, roughly. Timing on a 512MB 2.6G pentium, standard Redhat 9 gcc.

Simplifications possible..

AM I think you are missing the point:

I want to be able to use Tcl's support for parsing an arbitrary expression (that is why lexpr is not the solution for me)
I want to be able to avoid having to program the loops myself
I want to get the quickest response possible (without resorting to C or Fortran)

It is both an intellectual exercise and a matter of practical concern.

TV I think that I'm missing a profound point because there simply isn't an all to important point to this, unless I am led to belief what incredible range of interesting questions this can all lead to, and that that is the point.

My fortran knowledge is limited and from a trs-80 based (ancient) version and some university stuff (I recall [SPICE] was done in it), but the simple stuff I read in your example is done by my C program in about the same space, and more along the lines of

   command arguments

syntax tcl is famous for.

And without requiring special variable names (as fortran does) or special libraries (2 ines of macro could be put in a .h file easily).

And it works fine, as I showed even quite efficient. And at least insightfull.

Which is not something I can say about the code below, frankly I don't see any intellectual point to it except: stick with the subject you are trying to make something clear about. The only explicit subject I see raised in your introduction is a certain language construct to act on matrices (practically you deal only with vectors), allowing an operation to be automatically applied to all members or elements of the vector. That I understand.

Also, I understand you may want to execute a little tcl script on those members, which is without question a valid idea, and also I easily understand that in that case you are concerned with maybe not just allow expr to be used, maybe minding {} in the expr for efficiency, and preventing having to define procs (functions), maybe replacing them with a lambda equivalent.

I also understand easily you may want to define operations in a generic way for parts of a matrix (not vector, as you only use), which would be interesting challenge. Relative indices and boundary handling and all that. Like so many who claim that area because it is appealing to deal with it right, and powerfull, you don't even touch a solution for such things, nor did you directly claim to.

For numerical matrix operations, using a tcl expression is just as powerfull as a C expression, except of course in tcl you can change the expression run-time.

If you don't want to program the loop (I guess eastetics or readability issue), use a simple function with foreach arguments of something, I mean that's not very profound in tcl, unless you want to deal with the above.

When you mention speed and fortran yourself (a point I didn't make and understood perfectly well) I thought it was a decent, intellectually interesting and practical addition to the page idea to mention that people nowadays can often easily program in quite powerfull C, which can even be combined with tcl in runtime, on PC's and workstations, and that that gives incredible efficiency benefits. Three orders of magnitude, when I'm not mistaking, which in any measure is considerable. And it would not seem necessary to go steps back to such a relatively old fashioned language as fortran, when the same idea is, minus some syntactic sugar a full grown programmer probably wouldn't care much for, quite realizable and overseeable (in about the same code size), and for me excellent readability (like Pascal to me isn't much more readable then well-applied C).

When one wants to do it all in tcl, than the code below has relatively little to offer, and is an incredible hassle IMO.

I wasn't trying to invent a perpetuum mobile, and don't like it when I see a lot of claims made which strain my imagination as to what the real value is a but too far.

What you are without question at least alluding to is that you may want constructions/functions/tools to deal with arrays symbolically, which is of course interesting, but your example leaves me with nothing to realy see that happen, except something of the kind that fits in 10 decent lines of C code and repeated function/macro calling.

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

 }

AM An alternative implementation of "numset" is given below:

 #
 # Alternative implementation
 #
 proc ::numarray::numset2 { varname expression } {
    upvar $varname var

    #
    # Replace the substrings "$var" by a corresponding
    # list value - if the variable is a list
    #
    set _length_ 1
    set _count_  0
    set _vars_   {}
    set _empty_  {}

    foreach lv {lv0 lv1 lv2 lv3 lv4 lv5 lv6 lv7 lv8 lv9} {
       set $lv "_empty_"
    }

    while { [string first \$ $expression] >= 0 } {
       regexp {\$(\w+)} $expression ==> _v_
 
       #
       # If it is a list, then introduce a new variable
       # (Note: no more than 10 list variables right now)
       #
       upvar $_v_ $_v_
       if { [llength [set $_v_]] > 1 } {
          if { [lsearch $_vars_ $_v_] < 0 } {
             set replace "vv$_count_"
             unset lv$_count_
             upvar $_v_ lv$_count_
             lappend _vars_ $_v_
             incr _count_
          }
          set _length_ [llength [set lv$_count_]]
          #regsub {\$\w+} $expression "@$replace" expression
          set expression [string map [list \$$_v_ @$replace] $expression]
       } else {
          catch { upvar $_v_ $_v_ }
          #regsub {\$\w+} $expression "@$_v_" expression
          set expression [string map [list \$$_v_ @$_v_] $expression]
       }
    }

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

    #
    # Store the result in a separate variable - we should not
    # lose the old value yet!
    #
    set _result_ {}
    if { $_count_ == 0 } {
       lappend _result_ [expr $expression]
    } elseif { $_count_ == 1 } {
       foreach vv0 $lv0 {
          lappend _result_ [expr $expression]
       }
    } elseif { $_count_ == 2 } {
       foreach vv0 $lv0 vv1 $lv1 {
          lappend _result_ [expr $expression]
       }
    } elseif { $_count_ == 3 } {
       foreach vv0 $lv0 vv1 $lv1 vv2 $lv2 {
          lappend _result_ [expr $expression]
       }
    } else {
       foreach vv0 $lv0 vv1 $lv1 vv2 $lv2 \
               vv3 $lv3 vv4 $lv4 vv5 $lv5 \
               vv6 $lv6 vv7 $lv7 vv8 $lv8 \
               vv9 $lv9 {
          lappend _result_ [expr $expression]
       } 
    }
    set var $_result_
 }

As can be seen from the test results, the performance with this second version is not quite as good as the raw foreach loop. Miguel Sofer MS pointed out that the expression is not byte-compiled and suggested using a temporary proc to force it to be compiled. So this resulted in a third version:

 #
 # Third version 
 #
 proc ::numarray::numset3 { varname expression } {
    upvar $varname var

    #
    # Replace the substrings "$var" by a corresponding
    # list value - if the variable is a list
    #
    set _length_ 1
    set _count_  0
    set _vars_   {}
    set _empty_  {}

    foreach lv {lv0 lv1 lv2 lv3 lv4 lv5 lv6 lv7 lv8 lv9} {
       set $lv "_empty_"
    }

    while { [string first \$ $expression] >= 0 } {
       regexp {\$(\w+)} $expression ==> _v_
 
       #
       # If it is a list, then introduce a new variable
       # (Note: no more than 10 list variables right now)
       #
       upvar $_v_ $_v_
       if { [llength [set $_v_]] > 1 } {
          if { [lsearch $_vars_ $_v_] < 0 } {
             set replace "vv$_count_"
             unset lv$_count_
             upvar $_v_ lv$_count_
             lappend _vars_ $_v_
             incr _count_
          }
          set _length_ [llength [set lv$_count_]]
          #regsub {\$\w+} $expression "@$replace" expression
          set expression [string map [list \$$_v_ @$replace] $expression]
       } else {
          # 
          # Substitute the value
          #
          catch { upvar $_v_ $_v_ }
          #regsub {\$\w+} $expression "@$_v_" expression
          set expression [string map [list \$$_v_ [set $_v_]] $expression]
       }
    }

    set expression [string map {@ $} $expression]
 
    #
    # Create temporary procedures, so as to guarantee compilation 
    # to byte-code
    #
    if { $_count_ == 0 } {
       set var [expr $expression]
    } elseif { $_count_ == 1 } {
       proc ::numarray::lexpr {lv0} \
          "set _result_ {}
          foreach vv0 \$lv0 {
             lappend _result_ \[expr {$expression}\]
          }
          return \$_result_"
       set var [lexpr $lv0]
     } elseif { $_count_ == 2 } {
       proc ::numarray::lexpr {lv0 lv1} \
          "set _result_ {}
         foreach vv0 \$lv0 vv1 \$lv1 {
             lappend _result_ \[expr {$expression}\]
          }
          return \$_result_"
       set var [lexpr $lv0 $lv1]
    } elseif { $_count_ == 3 } {
       proc ::numarray::lexpr {lv0 lv1 lv2} \
          "set _result_ {}
          foreach vv0 \$lv0 vv1 \$lv1 vv2 \$lv2 {
             lappend _result_ \[expr {$expression}\]
          }
          return \$_result_"
       set var [lexpr $lv0 $lv1 $lv2]
     } else {
       proc ::numarray::lexpr {lv0 lv1 lv2} \
          "set _result_ {}
          foreach vv0 \$lv0 vv1 \$lv1 vv2 \$lv2 \
                  vv3 \$lv3 vv4 \$lv4 vv5 \$lv5 \
                  vv6 \$lv6 vv7 \$lv7 vv8 \$lv8 \
                  vv9 \$lv9 {
             lappend _result_ \[expr {$expression}\]
          }
          return \$_result_"
       set var [lexpr $lv0 $lv1 $lv2 $lv3 $lv4 $lv5 $lv6 $lv7 $lv8 $lv9]
    }
 }

AM With the following test code:

 #
 # Simple measurements
 #
 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 "numset:   [time {::numarray::numset a {$v+$w+$v*$w}} 10]"
 puts "testloop: [time {set a [testloop $v $w]} 10]"
 puts "numset2:  [time {::numarray::numset2 a {$v+$w+$v*$w}} 10]"
 puts "numset3:  [time {::numarray::numset3 a {$v+$w+$v*$w}} 10]"

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

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

I get these results:

 numset:   36055 microseconds per iteration
 testloop: 10859 microseconds per iteration
 numset2:  18213 microseconds per iteration
 numset3:  10293 microseconds per iteration

 numset:   3678 microseconds per iteration
 testloop: 1010 microseconds per iteration
 numset2:  1870 microseconds per iteration
 numset3:  1100 microseconds per iteration

 numset:   534 microseconds per iteration
 testloop: 108 microseconds per iteration
 numset2:  289 microseconds per iteration
 numset3:  234 microseconds per iteration

so effectively:

numset is four to five times slower than the tight loop
numset2 is two to three times slower (it gets better on long lists)
numset3 is just about as fast as the tight loop, for long lists and only two times as slow for short lists.

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.