Version 1 of Element-wise array operations

Updated 2004-10-21 08:43:24 by AM

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 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) an 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:

  • Arrays are considered to be 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

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