Arjen Markus (16 october 2020) In the Tclers' chat the other day, Steve Blinkhorn mentioned an article in Nature on the numPy package that is used for scientific and technical calculation in Python. The question he posed was whether it would be an idea to build an interface to numPy, so that we can use it from a Tcl program. One reason for this: numPy is a well-respected and well-known package and it would potentially make it easier for Tclers to get their work accepted.
Steve Landers told me of the tclpy package that enables you to call Python from within Tcl and vice versa (see Accessing Tcl and Python from one another for more information, as there are solutions too). So, this provided me a first step towards using numPy in a Tcl program. Well, "program" is a bit of a big word, as I have not tried anything serious yet. Still, it is a useful first step to investigate how this package can be used.
The preparations were simple: build the package and load it into Tcl. My Cygwin environment presented a few problems, though, as the makefile required a file python-config to provide the Python include files and libraries. After some experimenting I used the following command to build it:
gcc -o libtclpy0.3.so -shared generic/tclpy.c -I/usr/include/tcl8.6 -I/usr/include/python2.7 -ltcl8.6 -lpython2.7 -Wl,--export-all-symbols
and introduced a few missing macros in the source code:
#define PACKAGE_VERSION "0.3" #define PY_LIBFILE "libpython2.7.dll"
Then it was a matter of trying a few commands. Like from the documentation:
package require tclpy
py eval {import numpy as np}
py eval {def xx(x,y): return x+y}
puts [py call xx string1 string2]which diligently prints string1string2
Defining an array in the Python interpreter is a trifle more demanding:
set a [list 4 2 3] py eval "a = np.array(\[[join $a ,]\])"
to convert the data in the list into the right format for Python.
Now, however, the limitation: the tclpy package currently passes strings, not variables. This makes it a trifle difficult to manipulate variables on the Python whose names we pass from the Tcl side. For instance: to sort the numbers in the array we just created, this does not work:
puts [py call np.sort a]
Instead we need to convert the string "a" into the Python variable "a". Luckily Python does offer that possibility via the globals() and locals() function.
So what we can do:
py eval {def tolist(array): var = globals().get(array); return var.tolist()}
puts [py call tolist a]This works also with two-dimensional arrays
set a {[[4,2,3],[1,2,3]]}
py eval "z = np.array($a)"
puts [py call tolist z]prints:
{4 2 3} {1 2 3}Of course, this is not very efficient - with a large array or matrix you would need to convert everything to a long string, pass it to the Python interpreter, do the computations and pass the result back.
But the fact that it is so easy to do these first few steps makes it possible to think about the next steps: provide more direct transfer methods to pass large amounts of numerical data to and fro.
AM (2 november 2020) Here is a small program to test the performance of the interface. Note: this still uses the Python 2 version - in the meantime Aidan has updated the code to support Python 3, but I have no tested that yet.
# perform.tcl --
# Simple checks on the performance of the interface to Python
#
# Check: setting up matrices of different sizes
# Idea:
# Does it make a difference if you fill the matrix all at once or
# fill it one element at a time?
# What about filling it by row or by column?
# Matrix sizes: 10x10, 30x30 and 100x100
#
# Clear result:
# Setting the values in the whole matrix at once is 3x faster
# than setting the entries one by one
#
# Further checks:
# - Getting back the matrix
# - Solving a linear system (comparing with math::linearalgebra?)
#
proc setMatrix {size} {
set matrix {}
for {set r 0} {$r < $size} {incr r} {
set row {}
for {set c 0} {$c < $size} {incr c} {
lappend row [expr {rand()}]
}
lappend matrix "\[[join $row ,]\]"
}
set pythonMatrix "\[[join $matrix ,]\]"
py eval "matrix = np.array($pythonMatrix)"
}
proc setByElement {size} {
for {set r 0} {$r < $size} {incr r} {
for {set c 0} {$c < $size} {incr c} {
set rnd [expr {rand()}]
py eval "matrix\[$r,$c\] = $r"
}
}
}
proc getMatrix {size} {
py call tolist matrix
}
proc getByElement {size} {
for {set r 0} {$r < $size} {incr r} {
for {set c 0} {$c < $size} {incr c} {
py call tostr matrix $r $c
}
}
}
proc solveEquations {} {
py eval "x = np.linalg.solve(matrix,b)"
set x [py call tolist x]
}
set auto_path [concat . $auto_path]
package require tclpy
py eval {import numpy as np}
py eval {def tolist(array): aa = globals().get(array); return aa.tolist()}
py eval {def tostr(array,r,c): aa = globals().get(array); rr = int(r); cc = int(c); return aa[rr,cc]}
# Test:
# First set up the matrix on the Python side, then
# use the interface to fill the entries
# Use random numbers to make sure no string representation
# exists.
#
foreach size {10 30 100} counts {100 30 10} {
set size2 [expr {$size * $size}]
py eval "matrix = np.zeros($size2).reshape($size,$size)"
set result1 [time {setMatrix $size} $counts]
set result2 [time {setByElement $size} $counts]
#set result3 [time {setByRow $size} $counts]
#set result4 [time {setByColumn $size} $counts]
puts "Size: $size"
puts " Fill the matrix at once: $result1"
puts " Fill the matrix one entry at a time: $result2"
#puts " Fill the matrix one row at a time: $result3"
#puts " Fill the matrix one column at a time: $result4"
}
#
# Now retrieve the values of a matrix
#
foreach size {10 30 100} counts {100 30 10} {
set size2 [expr {$size * $size}]
py eval "matrix = np.random.random(($size,$size))"
set result1 [time {getMatrix $size} $counts]
set result2 [time {getByElement $size} $counts]
#set result3 [time {setByRow $size} $counts]
#set result4 [time {setByColumn $size} $counts]
puts "Size: $size"
puts " Retrieve the matrix at once: $result1"
puts " Retrieve the matrix one entry at a time: $result2"
#puts " Fill the matrix one row at a time: $result3"
#puts " Fill the matrix one column at a time: $result4"
}
#
# Solve a linear system
#
foreach size {10 30 100} counts {100 30 10} {
py eval "matrix = np.random.random(($size,$size))"
py eval "b = np.random.random($size)"
set result [time {solveEquations} $counts]
puts "Size: $size"
puts " Solve a linear system: $result"
}
#puts [solveEquations]
The output on my system, for easy reference:
Size: 10
Fill the matrix at once: 508.48 microseconds per iteration
Fill the matrix one entry at a time: 1149.97 microseconds per iteration
Size: 30
Fill the matrix at once: 2593.766666666667 microseconds per iteration
Fill the matrix one entry at a time: 10133.833333333334 microseconds per iteration
Size: 100
Fill the matrix at once: 33364.6 microseconds per iteration
Fill the matrix one entry at a time: 110351.6 microseconds per iteration
Size: 10
Retrieve the matrix at once: 76.88 microseconds per iteration
Retrieve the matrix one entry at a time: 337.51 microseconds per iteration
Size: 30
Retrieve the matrix at once: 731.1 microseconds per iteration
Retrieve the matrix one entry at a time: 3083.6 microseconds per iteration
Size: 100
Retrieve the matrix at once: 7080.7 microseconds per iteration
Retrieve the matrix one entry at a time: 35484.4 microseconds per iteration
Size: 10
Solve a linear system: 81.17 microseconds per iteration
Size: 30
Solve a linear system: 87.33333333333333 microseconds per iteration
Size: 100
Solve a linear system: 396.8 microseconds per iteration AM (4 november 2020) Here is a further experiment - solving a simple reaction-diffusion equation. The transport (diffusion) part is done on the Python side, but the reaction part is done via Tcl. This necessitates the transfer to and from Python of the whole solution, but it is the most variable part of such problems. So, a Tcl implementation is preferred.
Note:
During the development of the code below I was confronted with mysterious error messages from the Python parser/interpreter. These problems quite likely originated from my Tcl-based intuition on how the commands are interpreted and my poor understanding of Python. Still, the ability to express the finite difference formula for the transport step in such a compact form is encouraging.
# chkpde.tcl --
# Use NumPy and libtclpy to solve a simple PDE
#
# The PDE in question is a reaction-diffusion equation defined
# on a square grid. The equation:
#
# dC/dt = D nabla C - k C + 1
#
# The discretisation is straightforward:
#
# Cij(new) = Cij + Dt * (D/Dx^2 (Ci-1,j + Ci+1,j + Ci,j-1 + Ci,j+1 - 4 Cij) - k Cij + 1)
#
# Defined on a square grid and the boundary conditions on the four sides are
# simply: C = 0.
#
# D = 1.0e-3 m2/s
# Dx = 0.1/size m
# Dt = 0.1/size s
# k = 0.1 /s
# T = 10 s
#
# With the experience on the performance in mind, we send full matrices back and forth
# between Tcl and Python
#
proc setMatrix {name matrix} {
set content {}
foreach row $matrix {
lappend content "\[[join $row ,]\]"
}
set pythonMatrix "\[[join $content ,]\]"
py eval "$name\[:,:\] = $pythonMatrix"
}
proc setByElement {size} {
for {set r 0} {$r < $size} {incr r} {
for {set c 0} {$c < $size} {incr c} {
set rnd [expr {rand()}]
py eval "matrix\[$r,$c\] = $r"
}
}
}
proc getMatrix {name} {
py call tolist $name
}
proc calculateReaction {conc} {
global decay
set reaction {}
foreach row $conc {
set newrow {}
foreach col $row {
lappend newrow [expr {-$decay * $col + 1.0}]
}
lappend reaction $newrow
}
return $reaction
}
#
# Load the package
#
set auto_path [concat . $auto_path]
package require tclpy
py eval {import numpy as np}
py eval {def tolist(array): aa = globals().get(array); return aa.tolist()}
py eval {def tostr(array,r,c): aa = globals().get(array); rr = int(r); cc = int(c); return aa[rr,cc]}
#
# Try different grid sizes
#
set T 10.0
foreach size {5 10 30 100} {
puts "Size: $size"
set size2 [expr {$size * $size}]
set dt [expr {1.0 / $size}]
set dx [expr {1.0 / $size}]
set diff 1.0e-3
set decay 0.1
#
# Set the parameters
#
py eval "dt = $dt"
py eval "dx = $dx"
py eval "diff = $diff"
py eval "decay = $decay"
set rcentre [expr { $size / 2 }]
set ccentre [expr { $size / 2 }]
#
# Initial condition
#
py eval "matrix = np.zeros($size2).reshape($size,$size)"
py eval "dmatrix = np.zeros($size2).reshape($size,$size)"
py eval "reaction = np.zeros($size2).reshape($size,$size)"
#
# Loop over time
#
set t 0.0
while { $t < $T } {
#
# Set a time step
#
set reaction [calculateReaction [getMatrix "matrix"]]
setMatrix "reaction" $reaction
#py eval {dmatrix[1:-1,1:-1] = diff / (dx*dx) * (matrix[0:-2,1:-1] + matrix[2:,1:-1] + matrix[1:-1,0:-2] + matrix[1:-1,2:] - 4.0 * matrix[1:-1,1:-1])}
#py eval {dmatrix[1:-1,1:-1] = dmatrix[1:-1,1:-1] + reaction[1:-1,1:-1]}
py eval {dmatrix[1:-1,1:-1] = diff / (dx*dx) * ( matrix[0:-2,1:-1] + matrix[2:,1:-1] + matrix[1:-1,0:-2] + matrix[1:-1,2:] - 4.0 * matrix[1:-1,1:-1] ) + reaction[1:-1,1:-1]}
py eval {matrix[1:-1,1:-1] = matrix[1:-1,1:-1] + dt * dmatrix[1:-1,1:-1]}
set t [expr {$t + $dt}]
#
# Examine the midpoint
#
set centre [py call tostr matrix $rcentre $ccentre]
if { abs($t - int($t+0.5)) < 0.5 * $dt } {
puts "$t\t$centre"
}
#puts $t
}
}And the output from the above:
Size: 5 1.0 0.960599948 1.9999999999999998 1.8271874881419155 3.0000000000000004 2.6071575895203023 4.000000000000001 3.30769748349534 5.000000000000002 3.935698818258164 6.000000000000003 4.497697658329985 7.0000000000000036 4.999836700287916 8.000000000000004 5.447844946834249 9.0 5.847030867503363 9.999999999999996 6.202285768638312 Size: 10 0.9999999999999999 0.9561787941780185 2.0000000000000004 1.8209076395834303 3.0000000000000013 2.602822500607819 4.000000000000002 3.3096203144671 4.999999999999998 3.9481744379070376 5.999999999999995 4.524655800808933 6.999999999999991 5.044640537773083 7.999999999999988 5.513198185316464 8.999999999999984 5.934960816904726 9.99999999999998 6.31417603667851 Size: 30 0.9999999999999999 0.9531371154286179 2.0000000000000027 1.8154271946171865 2.999999999999999 2.5955291739697137 3.9999999999999956 3.301275936820813 4.999999999999992 3.9397486027667545 5.9999999999999885 4.51733657940186 6.999999999999985 5.039788170428689 7.999999999999981 5.512261607539829 8.999999999999979 5.939382226970672 9.999999999999975 6.325304268761792 Size: 100 1.0000000000000007 0.9520785288629056 2.0000000000000013 1.8135117052136498 2.99999999999998 2.5929296783229674 3.9999999999999587 3.2981409141647706 4.999999999999938 3.9362097981062054 5.9999999999999165 4.513523291350504 6.999999999999895 5.035844462033559 7.999999999999874 5.5083581162268365 8.999999999999853 5.93571804353677 9.999999999999831 6.322100126712278
arjen - 2020-11-06 08:04:42
Some additional information: Aidan has updated the code for Python 3 and Paul Obermeier has added support for CMake so that the library can be built for all manner of environments (notably the ones on Windows: Windows itself, Cygwin, MinGW)
PO - 2021-01-08
The tclpy package is now part of the BAWT 1.3.0 framework.