Profiling code in Smalltalk

The goal is to benchmark the accelerated Large Integer arithmetic, and track the costs by profiling. But let's first blog about the profiling facility itself.

This will be illustrated in Squeak, but applies to many Smalltalk dialects.

Squeak by default has a MessageTally class that is really easy to use, just evaluate code like this (my ongoing work on exponentiation):

MessageTally spyOn: [10000 timesRepeat: [103 raisedToInteger: 8000]].

How does it work? By setting a timer that will periodically interrupt the profiled Process, scanning the interrupted Process stack, and gather statistics about the calling sequence tree (it basically count the number of time that each method in the call tree has been interrupted). It's a gift from original Smalltalk-80, and is completely written at image side in pure Smalltalk, coexisting in the lively environment with the profiled code, and with other essential tools like Debugger. This is amazing enough and demonstrates the remarkable reflexion facilities of Smalltalk upon itself (introspection). I invite everyone to read that well written piece of code.

It also works remarkably well for analyzing pure Smalltalk code.

But for code which is dominated by long running primitives like Large Integer arithmetic, it does not work so well, because execution of primitives are atomic un-interruptible operations from the Smalltalk image point of view. The resut is that timer-based interruption is delayed and the profile cost attributed to the wrong message, the one following the primitive, and when profiling several processes, it's worse, see this mailing list message from Eliot Miranda, the main contributor of opensmalltalk-vm, or should I say the deus ex virtual machina.

Effectively, see the leaves at end of the report, it indicates that exponentiation time is dominated by addition! (no, we do not implement multiplication by repetitive additions).

 - 4631 tallies, 4644 msec.

**Tree**
--------------------------------
Process: (40) 58662: nil
--------------------------------
99.5% {4622ms} MessageTally class>>spyOn:reportOtherProcesses:
  99.5% {4622ms} MessageTally>>spyEvery:on:
    99.5% {4622ms} BlockClosure>>ensure:
      99.5% {4622ms} [] UndefinedObject>>DoIt
        99.5% {4622ms} SmallInteger(Integer)>>timesRepeat:
          99.5% {4620ms} [[]] UndefinedObject>>DoIt
            99.4% {4618ms} SmallInteger(Integer)>>raisedToInteger:
              99.4% {4618ms} SmallInteger(Integer)>>ternaryBinaryExponentationOf:
                96.8% {4495ms} LargePositiveInteger>>squared
                  |94.4% {4385ms} LargePositiveInteger>>squaredByFourth
                  |  |37.0% {1718ms} LargePositiveInteger>>*
                  |  |  |34.7% {1614ms} LargePositiveInteger(Integer)>>*
                  |  |  |  |34.7% {1614ms} LargePositiveInteger>>multiplyByInteger:
                  |  |  |  |  26.0% {1205ms} LargePositiveInteger>>digitMul22:
                  |  |  |  |    |17.3% {803ms} LargePositiveInteger>>+
                  |  |  |  |    |6.8% {316ms} LargePositiveInteger>>multiplyByInteger:
                  |  |  |  |  8.7% {405ms} LargeNegativeInteger(LargePositiveInteger)>>multiplyByInteger:
                  |  |  |  |    8.7% {405ms} LargeNegativeInteger(LargePositiveInteger)>>digitMul22:
                  |  |  |  |      5.7% {263ms} LargeNegativeInteger(LargePositiveInteger)>>+
                  |  |  |  |      2.2% {102ms} LargeNegativeInteger(LargePositiveInteger)>>multiplyByInteger:
                  |  |  |2.2% {104ms} primitives
                  |  |31.4% {1460ms} LargePositiveInteger>>squared
                  |  |  |24.9% {1155ms} LargePositiveInteger>>squaredByFourth
                  |  |  |  |9.1% {421ms} LargeNegativeInteger(Integer)>>bitShift:
                  |  |  |  |  |9.0% {419ms} primitives
                  |  |  |  |7.8% {363ms} LargePositiveInteger>>+
                  |  |  |  |4.5% {209ms} LargeNegativeInteger(LargePositiveInteger)>>*
                  |  |  |  |  |4.3% {201ms} primitives
                  |  |  |  |2.2% {100ms} LargePositiveInteger>>-
                  |  |  |  |  2.0% {94ms} primitives
                  |  |  |6.0% {278ms} LargePositiveInteger>>squaredByHalf
                  |  |  |  2.4% {110ms} LargePositiveInteger>>inplaceAddNonOverlapping:digitShiftBy:
                  |  |  |  1.7% {77ms} LargePositiveInteger>>squared
                  |  |  |    |1.6% {76ms} primitives
                  |  |  |  1.5% {69ms} LargePositiveInteger>>multiplyByInteger:
                  |  |14.6% {679ms} LargePositiveInteger(Integer)>>bitShift:
                  |  |8.6% {397ms} LargePositiveInteger>>+
                  |  |2.2% {104ms} LargePositiveInteger>>-
                  |2.2% {103ms} LargePositiveInteger>>squaredByHalf
                2.1% {98ms} primitives

**Leaves**
39.9% {1853ms} LargePositiveInteger>>+
23.6% {1098ms} LargeNegativeInteger(Integer)>>bitShift:
11.4% {529ms} LargePositiveInteger>>multiplyByInteger:
7.0% {325ms} LargeNegativeInteger(LargePositiveInteger)>>*
5.2% {243ms} LargePositiveInteger>>-
4.0% {186ms} LargePositiveInteger>>inplaceAddNonOverlapping:digitShiftBy:
3.6% {167ms} LargePositiveInteger>>squared
2.2% {102ms} SmallInteger(Integer)>>ternaryBinaryExponentationOf:

Since primitives are opaque at image side, it is necessary to introduce additional VM support for measuring such primitive cost. The simple idea for intercepting primitives is to instrument the VM at time of method activation. The timer-based interrupt is replaced by a polling of high resolution clock inside the VM. If the method is a primitive, and the sample interval is expired, then the currently active Process is recorded into profileProcess inst. var of the Interpreter. and the activated method (newMethod inst. var.) is recorded in the profileMethod inst. var., then the Semaphore of the profiler Process is signalled. This profileMethod and profileProcess can then be queried by the profiler thru primitives. The polling is also performed at each interrupt check - that is at each and every potential Process interruption point (method activation and long jump byte code). Since the active Process is recorded in profileProcess, before eventual Process switch, this avoid attribution of profiling cost to the newly running Process.

This support already exists and is covered by a set of four primitives accessible from the image:

• primitiveProfileSemaphore for setting the Semaphore to signal when a sample is ready (the profiler will wait on that semaphore).
• primitiveProfileStart for starting/stopping the profiling and setting the number of high resolution clock ticks between two sample
• primitiveProfilePrimitive for querying the profileMethod
• primitiveProfileSample for querying the profileProcess

The generated C code for these primitives and other support routines is available at opensmalltalk-vm, but it's better to browse the corresponding VMMaker Smalltalk source code from a Squeak image, because most of the VM is written in Smalltalk. For the curious, browse implementors of checkForInterrupts and internalPrimitiveResponse, or above mentionned inst. var. references.

The rest of profiling magic happens at image side, like with MessageTally. The profiler is not included in trunk images, but easily available at squeaksource3, or via SqueakMap.

It is named AndreasSystemProfiler in memory of Andreas Raab which was one of the major developer and inspiring leader of Squeak and one of the authors of this piece of code among many others. This time it's a gift from Ron Teitelbaum which agreed to release under MIT license in memory of Andreas (such re-license request was one of his suggestions).

The use of AndreasSystemProfiler is equally simple:

AndreasSystemProfiler spyOn: [10000 timesRepeat: [103 raisedToInteger: 8000]].

The report obtained is below:

Reporting - 45,077 tallies, 5,101 msec.

**Tree**
... snip ...
[                              99.94 (5,098)  UndefinedObject DoIt
[                                99.94 (5,098)  AndreasSystemProfiler class spyOn:
[                                  99.94 (5,098)  BlockClosure ensure:
[                                    99.94 (5,098)  [] AndreasSystemProfiler class spyOn:
[                                      99.94 (5,098)  AndreasSystemProfiler spyOn:
[                                        99.94 (5,098)  BlockClosure ensure:
[                                          99.94 (5,098)  [] UndefinedObject DoIt
[                                            99.94 (5,098)  SmallInteger(Integer) timesRepeat:
[                                              99.94 (5,098)  [[]] UndefinedObject DoIt
[                                                99.94 (5,098)  SmallInteger(Integer) raisedToInteger:
[                                                  99.94 (5,098)  SmallInteger(Integer) ternaryBinaryExponentationOf:
[                                                    96.92 (4,944)  LargePositiveInteger squared
[                                                      |94.24 (4,807)  LargePositiveInteger squaredByFourth
[                                                      |  |40.45 (2,063)  LargePositiveInteger *
[                                                      |  |  |37.78 (1,927)  LargePositiveInteger(Integer) *
[                                                      |  |  |  |37.78 (1,927)  LargePositiveInteger multiplyByInteger:
[                                                      |  |  |  |  29.27 (1,493)  LargePositiveInteger digitMul22:
[                                                      |  |  |  |    |15.47 (789)  LargePositiveInteger +
[                                                      |  |  |  |    |  |13.62 (695)  Integer digitMultiply:neg:
[                                                      |  |  |  |    |10.74 (548)  LargePositiveInteger multiplyByInteger:
[                                                      |  |  |  |    |  |9.12 (465)  Integer digitMultiply:neg:
[                                                      |  |  |  |    |1.11 (57)  LargePositiveInteger inplaceAddNonOverlapping:digitShiftBy:
[                                                      |  |  |  |  8.5 (434)  LargeNegativeInteger(LargePositiveInteger) multiplyByInteger:
[                                                      |  |  |  |    8.5 (434)  LargeNegativeInteger(LargePositiveInteger) digitMul22:
[                                                      |  |  |  |      4.16 (212)  LargeNegativeInteger(LargePositiveInteger) multiplyByInteger:
[                                                      |  |  |  |        |4.11 (210)  Integer digitMultiply:neg:
[                                                      |  |  |  |      4.08 (208)  LargeNegativeInteger(LargePositiveInteger) +
[                                                      |  |  |  |        3.9 (199)  Integer digitMultiply:neg:
[                                                      |  |  |1.12 (57)  Integer digitMultiply:neg:
[                                                      |  |26.74 (1,364)  LargePositiveInteger squared
[                                                      |  |  |19.34 (987)  LargePositiveInteger squaredByFourth
[                                                      |  |  |  |6.94 (354)  LargePositiveInteger +
[                                                      |  |  |  |  |2.77 (142)  Integer digitMultiply:neg:
[                                                      |  |  |  |  |1.41 (72)  Integer digitAdd:
[                                                      |  |  |  |4.78 (244)  LargeNegativeInteger(LargePositiveInteger) *
[                                                      |  |  |  |  |4.01 (205)  Integer digitMultiply:neg:
[                                                      |  |  |  |4.76 (243)  LargePositiveInteger(Integer) bitShift:
[                                                      |  |  |  |  |2.46 (125)  Integer digitMultiply:neg:
[                                                      |  |  |  |  |1.32 (68)  Integer digitSubtract:
[                                                      |  |  |  |1.27 (65)  LargePositiveInteger -
[                                                      |  |  |  |1.26 (64)  LargePositiveInteger squared
[                                                      |  |  |6.99 (357)  LargePositiveInteger squaredByHalf
[                                                      |  |  |  3.68 (188)  LargePositiveInteger inplaceAddNonOverlapping:digitShiftBy:
[                                                      |  |  |    |3.42 (175)  Integer digitMultiply:neg:
[                                                      |  |  |  1.82 (93)  LargePositiveInteger squared
[                                                      |  |  |    |1.59 (81)  Integer digitMultiply:neg:
[                                                      |  |  |  1.2 (61)  LargePositiveInteger multiplyByInteger:
[                                                      |  |  |    1.14 (58)  Integer digitMultiply:neg:
[                                                      |  |14.42 (736)  LargePositiveInteger(Integer) bitShift:
[                                                      |  |  |10.68 (545)  Integer digitMultiply:neg:
[                                                      |  |  |2.36 (121)  Integer digitAdd:
[                                                      |  |7.4 (377)  LargePositiveInteger +
[                                                      |  |  |2.51 (128)  Integer digitMultiply:neg:
[                                                      |  |  |2.14 (109)  Integer bitShiftMagnitude:
[                                                      |  |3.9 (199)  LargePositiveInteger -
[                                                      |  |  1.67 (85)  Integer digitSubtract:
[                                                      |2.66 (136)  LargePositiveInteger squaredByHalf
[                                                      |  2.19 (112)  LargePositiveInteger multiplyByInteger:
[                                                      |    2.18 (111)  Integer digitMultiply:neg:
[                                                    2.43 (124)  Integer digitAdd:

**Leaves**
63.12 (3,220)  Integer digitMultiply:neg:
8.28 (423)  Integer digitAdd:
6.68 (341)  Integer digitSubtract:
6.61 (337)  Integer bitShiftMagnitude:


The tree is more detailed because sample interval has been reduced to about 100 µs, which is not possible with image-side timer-based polling. We don't get exactly 10 samples per millisecond (45 000 tallies for 5 000 ms), because long primitives execution time may cause some pause in polling, and drift the sampling interval I presume. Fortunately, divide and conquer split the time spent in primitives in tiny fragments, which is good for this profiling strategy.
The tree is still not perfect because the primitives sometimes get attached to the wrong caller (like + sending digitMultiply:neg: ???).
But the leaves are much better: we see that the VM is effectively spending most time in multiplication, even if time spent in addition/subtraction/shifting is not neglectable. That's more realistic!
 

Accelerating Large Integer arithmetic in Squeak Smalltalk

I recently introduced some changes in Kernel package of Squeak development trunk http://source.squeak.org/trunk.html so as to accelerate Large Integer arithmetic involving huge numbers.

The operations boosted by this set of changes are:

• multiplications *
• divisions / // \\ quo: and rem:
• squaring squared
• square root sqrt sqrtFloor

All these changes can easily be made optionals. There is only one mandatory change in the Kernel package, the introduction of an indirection multiplyByInteger: and divideByInteger: in Integer so that subclass can define their own hook. Default behaviour of the hook is to call the schoolbook primitives primDigitMultiplyNegative via digitMultiply:neg: and primDigitDivNegative via digitDiv:neg:. I believe that Pharo and Cuis Smalltalk should adopt these basic changes, as the penalty is really neglectable, and make the accelerated arithmetic a package on its own. I've also inquired about Visualworks, but the changes required to the base classes are a bit more invasive.

For accelerating those operations, I have used these classical divide and conquer algorithms:

• Karatsuba multiplication in digitMul22:
• Toom-Cook multiplication in digitMul23: and digitMul33:
• the Fast Recursive Division - MPI-I-98-1-22 from Christoph Burnikel & Joachim Ziegler in digitDiv21: digitDiv32: and driver digitDivSplit:
•  and beautiful Karatsuba Square Root from Paul-Zimmerman - Rapport de recherche n° 3805
  
  lk*jFiGhFg6=f42/eHeW>4@/eH+4YW􏰓4dc`IRWGCba2aFJ`
  

Here, the 22, 23, 33, 21, and 32 in selector names tells in how many parts the operands (receiver and argument) are divided, 22 means two parts each (Karatsuba).

The multiplyByInteger: method is refined in LargePositiveInteger and serves as a dispatcher to the available variants, based on heuristics on the operands digitLengh. Remind that those digits are 1-byte long (8 bits) in traditional Smalltalk, though the LargeIntegersPlugin of opensmalltalk-vm uses 32-bits digits under the hood - that is why the split is made at four bytes boundary (notice the bitClear: 2r11 used in most digit-splitting methods).

When the operands sizes are unbalanced (the larger is twice longer than the smaller or more), a digitMulSplit: is used with a strategy of splitting the larger in chunks of about 1.5 times the smaller. I have tried other combinations, 1-1 and 1-2, and measured that 1-1.5 was faster, thus the addition of a digitMul23: variant of Toom-Cook. It's interesting to notice how I reconstruct the result from collected parts in O(N) by in-place modification. Interlacing is used so as to be sure to get non-overlapping digits (see inplaceAddNonOverlapping:digitShiftBy:). My first iteration did use a O(N log N) of the form (a3 * x + a2)*x^2 + (a1 * x + a0), where multiplication by x correspond to bitShift: of the appropriate digit length, but each operation means allocating memory, moving bytes, and ends up in non neglectable cost noticeable in profiling...

The 3-way Toom-Cook implemented in digitMul33: is a variant from Marco Bodrato & Alberto Zanoni, and it is interesting to read What About Toom-Cook Matrices Optimality? on this subject.

Squaring also use Karatsuba and Toom-Cook, and these are asymetrical versions which outperform the symetrical ones. Specifically, squaredByFourth use the 4-way Toom-Cook variant of  Jaewook Chung and M. Anwar Hasan - Asymmetric Squaring Formulae   https://www.lirmm.fr/arith18/papers/Chung-Squaring.pdf.

The new square root algorithms introduces a new message sqrtRem which answers an Array with truncated square root and remainder for almost the same price
self sqrtRem first squared + self sqrtRem last = self
 
This is just a starting activity, there is more to do: I have played with exponentation (raisedToInteger:), cubing (cubed) and cube root (cbrtRem as a generalization of sqrtRem) and maybe I'll blog about them, or publish if worth. There is also a need to auto-adapt the thresholds hardcoded here and there, but this is quite involved to get it right, benchmarking is an art.

Why bother with those classical algorithms without much innovating, when the state of the art GNU Multi Precision library (GMP) is far more advanced, optimized, fast and ready to use? Because I have an interest in supporting ArbitraryPrecisionFloat, and more than all, because it's very easy and fun to experiment in Squeak. It's fascinating how straight forward it is to transcribe those algorithms in Smalltalk and I believe that researchers should prototype new algorithms in such a language before paying the ultimate tribute to optimization. If you find that extending GMP is equally fun, which means programming in C, dealing with memory optimization, edit-compile-run-segfault cycles, then OK, you're even more a weirdo than I am!

 

Swimming against the tide

This is an old article of 2018 that I forgot to publish, about how argh-woops-dang the programming experience can be sometimes, maybe the title should have been swimming in bitumen.

HDF5 is an essential format for exchanging large chunks of data. I don't much like it. It's overly complex. But not having an interface to it is a severe handicap for a language aiming to be generalist, and crunching data is rather a generalist task these days.

Of course, I could be using Python like anyone else on earth, but I would prefer some more productive and lively environment. So a few months ago, I've started to write such an interface for exchanging data between HDF5 and Smalltalk (Visualworks for a beginning). I'll publish on public store when having minimal usable core features, and it will be licensed MIT. But it's not yet ready. It's an activity taken on my free time.

I've started development in windows because having a not-supported-by-apple library working on Macosx is rarely an "out of the box" experience. This week end, I decided to give it a kick and to switch back to my preferred platform.

Building the library with the cmake build configuration provided by the hdfgroup was fairly easy. Just follow the instructions from http://support.hdfgroup.org/HDF5/release/cmakebuild.html.... and produce what you don't need: static libraries.

Building the dynamic libraries has been more involved, because neither setting the STATIC_ONLY=NO nor NO_MAC_FORTRAN=YES in build-unix.sh did the trick, whatever the recommendations found in http://support.hdfgroup.org/HDF5/release/chgcmkbuild.html.
The worse is that I don't know why. Did I forget to remove the previous build?
It's certainly not the right way, but brute force hardcoding of set(ADD_BUILD_OPTIONS "${ADD_BUILD_OPTIONS} -DBUILD_SHARED_LIBS:BOOL=ON") in the various cmake files produced the libhdf5.dylib that I wanted...

Or not... These were x86_64 libraries.  With my dated Smalltalk environment, I wanted an i386, or better a fat 96bits universal version. Not a problem. As usual, Stackoverflow has THE solution: http://stackoverflow.com/questions/5334095/cmake-multiarchitecture-compilation. I stupidely passed option "-DCMAKE_OSX_ARCHITECTURES=x86_64;i386" to ctest, which was not the right incantation, and without any flashing illumination, but rather after laborious attempts, finally added it as build option inside CTestScript.cmake

set (BUILD_OPTIONS "${BUILD_OPTIONS} \"-DCMAKE_OSX_ARCHITECTURES=x86_64;i386\" "). 

This was the right syntax, but not the right way...
My move was anticipated for a long time by another chess player, see those lines in
 hdf5-1.10.1/config/cmake_ext_mod/ConfigureChecks.cmake:

if (APPLE)
  list (LENGTH CMAKE_OSX_ARCHITECTURES ARCH_LENGTH)
  if (ARCH_LENGTH GREATER 1)
    set (CMAKE_OSX_ARCHITECTURES "" CACHE STRING "" FORCE)
    message(FATAL_ERROR "Building Universal Binaries on OS X is NOT supported by the HDF5 project. This is"
    "due to technical reasons. The best approach would be build each architecture in separate directories"
    "and use the 'lipo' tool to combine them into a single executable or library. The 'CMAKE_OSX_ARCHITECTURES'"
    "variable has been set to a blank value which will build the default architecture for this system.")
  endif ()
  set (${HDF_PREFIX}_AC_APPLE_UNIVERSAL_BUILD 0)
endif ()

Ouch! I don't care much of the best approach. I just want something that works! 
And I don't even agree. The best approach is to follow the guidelines of the target platform. 
Everything else is swimming against the tide...
For the time being, it's me who is swimming against the tide, too far away from the mainstream.
Even if I don't build for universal but just i386 architecture, the option will be blanked.
So what are these technical reasons exactly? Where can I read about that?

Without a clue, I started to hack those lines if (ARCH_LENGTH GREATER 1000), and also  set (${HDF_PREFIX}_AC_APPLE_UNIVERSAL_BUILD 1) but of course, as advertised this did not work. OK, with -VV verbose verbose option, we know that the link failed because ZLib and SZip were compiled for x86_64 only too. Untar-ing the provided ZLib.tar.gz and SZip.tar.gz, modifying the ConfigureChecks.cmake inside which had the very same guard as above, re-taring-re-gzipping ... equally failed! Same problem, the external libraries were NOT compiled with universal support. Adding the DCMAKE_OSX_ARCHITECTURES=x86_64;i386 to various CMAKE_ARGS definitions that I found did not help either...

Hacking randomly like this is vain and probably doomed to fail, but how are you supposed to understand what happens and where it comes from exactly among those thousand lines scattered all around the .cmake files?

You think I'm exagerating? 

find . -name '*cmake*' | wc -l
-> 101

wc -l `find . -name '*cmake*'`
-> 23729

More than twenty thousand lines in a hundred files! With all that stuff, HDF5 should certainly compile on your android tablet and your connected fridge! But not on my Mac...

If you count the generated build directory, this roughly doubles.
But if you remove the duplicates (and blanks) it's only half:

cat `find . -name '*cmake*'` | sort |  uniq | wc -l
-> 12799

Believe me now?

Back to the technical reasons they spoke about, I now have a clue: could it be the choice of the right tools? Ugly scripts, ugly environment variables, ugly macros, duplicated code and no debugger for getting a chance to dive inside the machinery. Welcome in the 21st century.

That's what sucks with configure/cmake and all those meta levels that generates configuration files for yet another tool (make): you gradually loose control on the lower levels where things has to be solved. The problem is fairly simple: compile/link ZLib and SZip with the universal build flags. The solution of going across all these layers is really involved, even if I don't even want to make it work for others, but just for me!

I could have talked about Smalltalk and the implementation, but half the week-end is gone in smoke now without the slightest progress.

Slow HDF5 progress

I've been struggling with DLLCC, and found another bizarre behavior that I reported on S.O., http://stackoverflow.com/questions/49564024/isnt-pointer-type-checking-disabled-in-dll-c-connect-and-is-that-ok. I then took some time to reread the F... manual.

I think that I have deeper understanding of DLLCC now and published a new patch using double dispatching on the public store. This enables passing a Smallapack.CArrayAccessor as argument to the external function call. What's this thing (the name is not so nice)? It's a proxy to some C data (on heap), somehow like a CPointer (a CType plus a pointer on C data), but with two differences: 

1. it carries a length, and thus enables safe access of contents from within the Smalltalk world. If you think of it two minutes, every pointer should be bounded, that's clearly the way to go. Alas, in 2018 we still are at assembler level with those C API, so the pointer coming from external world still are unbounded (a pity for safety...).
2. it is using 1-based index and enables handling from within Smalltalk world like any other collection without too much mindstorm.

With the time spent on DLLCC, the progress was slow on HDF5, I just corrected the H5TString transfer. The documentation was not clear about fixed size null terminated strings. I first thought (or read?) that the terminating null was mandatory. But then the type H5T_C_S1 - a C type string of fixed size 1 - makes no sense! So I speculated that an extra byte was allocated for the terminator... This was wrong! If not enough place is allocated, then the null terminator is omitted. The difference of NullTerminated and NullPadded HDF5 Strings is then germane: the former ends before first null, the later at last non null, but in neither case null is mandatory...

Along with the CArrayAccessor change, I can now explore the contents of a HDF5 file created by Matlab save -v7.3. a Small step in the right direction.

HDF5 interface requires a patch in DLLCC

Apologies to those who would have tried the HDF5 interface (if any): I completely forgot that I had to patch argument coercion in DLL/C-Connect to let it work.

I suspect that this is a bug of DLLCC, see http://stackoverflow.com/questions/49544642/why-cant-i-pass-an-uninterpretedbytes-to-a-void-thru-dll-c-connect

I have published a new bundle with the override. Alternate solution would be to either pretend that the function can take an _oopref* buffer argument, or to go thru copyToHeap copyFromHeap: complications that I absolutely do not want (HDF5 data can be big, thus I do not wish any extra/unecessary copy).

Interfacing Smalltalk with HDF5

I am going to speak about connecting Smalltalk with external (scientific) data again.

This time, it's not with Matlab, but HDF5, the hierarchical data format supported by http://www.hdfgroup.org 

Every language targetting science/engineering niche must have an interface to HDF5. That is the case of Python and Matlab to only cite two. And you know that I'd like to promote the usage of Smalltalk in this area too. So lets do it: here comes the HDF5 bundle in Cincom public store.

If interfacing with Matlab mat-file format was a piece of cake, HDF5 is much more involved. First, because HDF5 is like a file system. There are recursive named Group of data, like directories. And groups are not necessarily arranged as a tree, but can form arbitrary graphs (circular) thanks to links. Second, because HDF5 comes with a type system. It can hold arbitrary types, whether atomic (integer, floating point,bit fields) or composite (structure, arrays). The types can even be references to other objects, which means that it is sufficiently general to describe heterogeneous collections of dynamically typed data.

There are two kind of way to store data: the first is Dataset which are named entries in Groups (a bit like files in directories). a Dataset is a rectilinear multiple dimension array (like the MultipleDimensionArray that I recently promoted in Visualworks public store and Squeak STEM). The second is Attributes. All named entries (Group, Dataset and named types) can have attributes, which are kind of string-key arbitrary-value pairs. Attributes generally hold meta data. They lack the ability to perform the read/write of sub-regions of data. We somehow can compare attributes to the property list of Morphic.

For additional complexity from the interfacing point of view, arrays can be of variable length. Thus the buffer for holding complex data cannot be preallocated from within Smalltalk before the data transfer, but has to be allocated on the fly by the HDF5 library. Thats potentially means either memory leaks or dangling references to freed memory, or the two if the programmer worked too late at night.

HDF5 comes with lot of documentation including reference, user guide, tutorial and examples if you want to learn more.

Among the many features of HDF5, let's focus of some of the most essential:

• the ability to perform type transformations during transfer operations. This enables language interoperability, since the type system is flexible enough to accomodate many languages (if not all thru C interface).
• the ability to read/write sub-regions of dataset. This enables handling of huge data, the whole dataset does not have to reside in memory.
• it scales well (or at least it can if used adequately).

That being said, this comes with a price: HDF5 is complex, and like the manipulated data, the API is bigger than big too.

It's clear that the target languages that HDF5 creators had in mind were statically typed. So the mapping to dynamically types objects is not straightforward, nor optimized, especially for the compound type (structure). The transfer necessarily involves two steps: the transfer or raw data, followed by a transformation to Smalltalk data for read, et vice et versa for write. For handling all cases, including references, a visitor pattern will be necessary (references can be cyclic too).

For huge data, the way to go is to use proxy to HDF5 objects. That means giving minimal behavior to the proxys, at least for storing/retrieving whole data or subregions. Application specific behavior should better lie in application specific objects, and those objects would use a generic HDF5 proxy. Since data structure is hierarchical (recursive), an application specific visitor should construct the application specific object graph. The early implementation that I just published is not yet there. It does not even use a visitor, but hardcode an arbitrary visit in the HDF5 proxies. It's more a minimal proof of concept at this stage, but a usable proof of concept though.

So, how do we use it? For creating a HDF5 file, first do this:

out:= H5File create: 'foo.h5'.
H5Dataset createPath: 'float1' parent: out rootGroup value: 1.3e0.
H5Dataset createPath: 'double2' parent: out rootGroup value: Double pi.
H5Dataset createPath: 'int3' parent: out rootGroup value: -357.
out close.

The close operation is not strictly necessary, because I implemented a registry of opened hdf5 objects with auto-close facility when the entries are reclaimed. But forcing a close fushes the file, otherwise the close will be delayed until all opened HDF5 entities have been reclaimed by garbage collector.

For reading, it's like this:

in:= H5File readOnly: 'foo.h5'.
float1 := (in / 'float1') value.
double2 := (in / 'double2') value.
int3 := (in / 'int3') value.
in close.

This example is a bit poor. We can also create/query groups, attributes and handle more complex data, thanks to MultipleDimensionArray, RawArray for atomic data types, the CArrayAccessor of Smallapack for arbitrary data. 

For writing Smalltalk objects on HDF5 files, one must define these 3 essential selectors:

• h5mem returns a buffer containing the data to transfer and usable by HDF5 API (must be compatible with DLL-C-Connect interface)
• h5type gives a HDF5 type description of buffer contents
• h5space gives  a HDF5 description of buffer layout (dimensions of the dataset)

Last note: I've used HDF5 version 1.8.10 for the interface. Unfortunately, HDF5 use macros for enabling backward compatibility, but the way DLL/C-connect works, we will have to subclass the H5Interface in order to support various versions.

That's all for this post, I may add more implementation details later on.

MatFileReader is ported to Squeak

I've ported the MatFileReader package initially written in Visualworks to Squeak. It should work unchanged in Pharo but I did not test. Usage in Squeak is simply:

    (MatFileReader on: (StandardFileStream readOnlyFileNamed: 'test.mat') binary) decode.

This will answer a Dictionary with workspace variable names as keys and MxArray as values which are essentially stubs by now. With a visitor pattern or by subclassing MatFileReader, it is possible to map the values to other Smalltalk classes.

Code repository is on squeaksource, http://www.squeaksource.com/STEM.html

    MCHttpRepository
        location: 'http://www.squeaksource.com/STEM'
        user: ''
        password: ''.

STEM stands for Smalltalk Tools for Engineering and Mathematics (or Squeak Tools if you're tainted). The name is short, close to the main acronym for Science Technology Engineering Mathematics, and otherwise have quite many meanings (http://en.wikipedia.org/wiki/Stem) . I think I will use it by now.

I've not assembled any MonticelloConfigurationMap nor prepared any Metacello ConfigurationOfSTEM, so by now you'll have to install all the packages in this repository by yourself.

Note that STEM is not a concurrent of Polymath (http://github.com/PolyMathOrg/PolyMath). Polymath is for Pharo, and I don't yet use Pharo, I'm far more fluent in Squeak. License is MIT, and volunteers for porting/integrating into Polymath are welcome if the library sound interesting enough. I just ask to have a little respect for authorship, if ever the environment still permits it.