Element-wise array operations

Arjen Markus (21 october 2004) I was facing the problem of having to read a file with a lot of floating-point data. The file contained a large set of timeseries (organised as data per time) and I needed to determine the average of each timeseries. So, in a similar fashion as Numerical array operations, I thought it should be possible to do this in a "generic" way:

   Loop over all the records of the file:
      Sum the data per time
   Divide all the sums by the number of times
   Output

Now, in Playing APL RS showed how to do this with (in principle) a sequence or array of arbitrary dimensions (in his case: using newlines to separate the levels). I wanted a bit more flexibility and performance though.

Here is my experiment:

The sequences are implemented as nested lists
If a list contains only one element, this is the end of the recursion
To avoid lots of very similar code, I use a template to generate specific procedures

Missing from this:

Unary operations
Typical operations from linear algebra, such as matrix multiplication

For instance:

{1 2 3 4 5 6} is a one-dimensional sequence (a one-dimensional array in C/Fortran)
{ {1 2} {3 4} {5 6}} is a two-dimensional sequence
{ { {1 2} {3 4} } { {5 6} {7 8} } is a three-dimensional sequence
and so on

The binary operations just combine the individual elements:

   {1 2 3 4 5 6} + {1 2 3 4 5 6} ==> {2 4 6 8 10 12}

 # arrayop.tcl --
 #    Element-wise array operations (arrays are implemented as
 #    lists or lists of lists, to any depth needed)
 #
 #    Note:
 #    Because the various operations require very similar code
 #    to handle arbitrary dimensions, I use a template to generate
 #    the actual code "on the fly"
 #

 namespace eval ::Arrays {

     # Create the namespace

 set CODE {
 # Lscalar_<op> --
 #     Perform a binary operation on a scalar and an array
 # Arguments:
 #     op1         Scalar in question
 #     op2         Array (nested list)
 # Result:
 #     Array with updated values
 #
 proc ::Arrays::Lscalar_<op> {op1 op2} {
     set v $op1
     if { [llength [lindex $op2 0]] <= 1 } {
         set result {}
         foreach w $op2 {
             lappend result [<expr>]
         }
     } else {
         set result {}
         foreach w $op2 {
             lappend result [Lscalar_<op> $v $w]
         }
     }
     return $result
 }
 # Rscalar_<op> --
 #     Perform a binary operation on an array and a scalar
 # Arguments:
 #     op1         Array (nested list)
 #     op2         Scalar in question
 # Result:
 #     Array with updated values
 #
 proc ::Arrays::Rscalar_<op> {op1 op2} {
     set w $op2
     if { [llength [lindex $op1 0]] <= 1 } {
         set result {}
         foreach v $op1 {
             lappend result [<expr>]
         }
     } else {
         set result {}
         foreach v $op1 {
             lappend result [Rscalar_<op> $v $w]
         }
     }
     return $result
 }
 # Array_<op> --
 #     Perform a binary operation on two arrays
 # Arguments:
 #     op1         First array (nested list)
 #     op2         Second array (nested list)
 # Result:
 #     Array with updated values
 #
 proc ::Arrays::Array_<op> {op1 op2} {
     if { [llength [lindex $op1 0]] <= 1 } {
         set result {}
         foreach v $op1 w $op2 {
             lappend result [<expr>]
         }
     } else {
         set result {}
         foreach v $op1 w $op2 {
             lappend result [Array_<op> $v $w]
         }
     }
     return $result
 }
 # <op> --
 #     Perform a binary operation on two arrays or scalars
 # Arguments:
 #     op1         First array (nested list) or scalar
 #     op2         Second array (nested list) or scalar
 # Result:
 #     Array with updated values
 #
 proc ::Arrays::<op> {op1 op2} {
     if { [llength $op1] <= 1 } {
         return [::Arrays::Lscalar_<op> $op1 $op2]
     }
     if { [llength $op2] <= 1 } {
         return [::Arrays::Rscalar_<op> $op1 $op2]
     } else {
         return [::Arrays::Array_<op> $op1 $op2]
     }
 }
 #
 # End of template
 }

 # add, sub, mult, div, max, min --
 #     Binary operations on arrays
 #
 eval [string map {<op> add  <expr> "expr {$v+$w}"} $CODE]
 eval [string map {<op> sub  <expr> "expr {$v-$w}"} $CODE]
 eval [string map {<op> mult <expr> "expr {$v*$w}"} $CODE]
 eval [string map {<op> div  <expr> "expr {$v/$w}"} $CODE]
 eval [string map {<op> max  <expr> "expr {$v>$w?$v:$w}"} $CODE]
 eval [string map {<op> min  <expr> "expr {$v<$w?$v:$w}"} $CODE]

 } ;# End of namespace

 # mkArray --
 #     Make an array with a uniform value
 # Arguments:
 #     ndim        "Toplevel" dimension
 #     args        Inferior dimensions and the value to be assigned
 # Result:
 #     Array with updated values
 #
 proc ::Arrays::mkArray {ndim args} {
     puts "mkArray: $ndim -- $args"
     if { [llength $args] < 1 } {
         error "At least two arguments required - dimension and value"
     }
     if { [llength $args] == 1 } {
         set v [lindex $args 0]
     } else {
         set v [eval mkArray $args]
     }

     set result {}
     set i      $ndim
     while { $i > 0 } {
         lappend result $v
         incr i -1
     }

     return $result
 }

 # test --
 #
 # Test the array construction
 #
 console show
 puts "2x3 array of 1's: [::Arrays::mkArray 2 3 1]"
 puts "5x3x2 array of x's: [::Arrays::mkArray 5 3 2 x]"

 set array1 { {1 2 3} {4 5 6} }
 set array2 { {6 5 4} {3 2 1} }
 set array3 { { {6 5 4} {3 2 1} } { {6 5 4} {3 2 1} } }
 set array4 { { {1 2 3} {4 5 6} } { {1 2 3} {4 5 6} } }

 puts "Sum: $array1+$array2 ="
 puts "     [::Arrays::add $array1 $array2]"
 puts "Sum: $array3+$array4 ="
 puts "     [::Arrays::add $array3 $array4]"
 puts "Min: min($array3,$array4) ="
 puts "     [::Arrays::min $array3 $array4]"
 puts "Min: min($array3,4) ="
 puts "     [::Arrays::min $array3 4]"

Category Numerical Analysis