## 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 ="