Liam Healy

05:19 pm - GSLL using FSBV

The GNU Scientific Library for Lisp (GSLL) now makes use of Foreign Structures By Value (FSBV) to define functions that pass or return complex scalars, which are represented as structures in GSL. The tests in gsll-tests will test some of these functions, and they all pass.

Since my last posting, there have been some improvements to FSBV; among them, the copying between Lisp and foreign code happens automatically in the definition with either fsbv:defcfun or fsbv:foreign-funcall. Since there is some preparation involving structures needed for libffi, and this need be done only once per function when it is defined, fsbv:defcfun makes a closure separately when it is loaded and then the invocation of the function just calls that closure on the specific arguments. If fsbv:foreign-funcall is called with the symbol used for the Lisp name of the fsbv:defcfun, it will also use that closure. GSLL takes advantage of that efficiency. Since I think FSBV won't change too much now (at least in the interface) I have tagged it with "0.1" release.

For those GSLL users who have no need for complex functions and don't want to bother installing libffi and FSBV, GSLL will compile and load fine without them; of course, those functions won't be available.

09:23 pm - Foreign structures by value

I have written Foreign Structures By Value which will allow calling foreign functions that either take or return structures by value instead of by reference (i.e., with a pointer). It requires CFFI and libffi, and should work in any environment where those are both supported.

For example, the complex type is defined in the GNU Scientific Library (GSL) as

typedef struct
  {
    double dat[2];
  }
gsl_complex;

The function gsl_complex_conjugate takes and returns a structure of this type by value,

gsl_complex gsl_complex_conjugate (gsl_complex z);

If we define

(fsbv:defcstruct (complex :constructor complex :deconstructor (realpart imagpart))
  (dat :double :count 2))
(defun complex-conjugate (complex-number)
  (fsbv:with-foreign-objects ((gslin 'complex complex-number))
    (object
     (fsbv:foreign-funcall "gsl_complex_conjugate" complex gslin complex)
     'complex)))

(complex-conjugate #c(3.0d0 4.0d0))
#C(3.0 -4.0)

This is making its way into GSLL but in the meantime should be usable for other projects.

09:00 pm - Regression (unit) test in Lisp, including float tests

In the course of writing GSLL, I found that I needed some software for regression (unit) tests. I mentioned before that I had made use of lisp-unit, augmenting with a test of floating point equivalence. Tom Hermann has similar needs; he saw what I had done and merged it with his own, better, floating point test. He has put this together into an augmented lisp-unit. The git repository has been set up with the permission of Christopher Riesbeck, the author of lisp-unit. Instead of including my extended lisp-unit now, GSLL (specifically, gsll-tests) makes use of this separate repository. Other projects could similarly benefit.

06:15 pm - Passing a complex number to a foreign function

GSL has a number of functions that take complex scalars as arguments. It defines them with either

typedef struct
  {
    double dat[2];
  }
gsl_complex;

typedef struct
  {
    float dat[2];
  }
gsl_complex_float;

I thought maybe there was a way around this. Could it be that the real and imaginary parts are simply passed sequentially as if they were two successive arguments to the function? In fact, this does work for numbers of type (complex double-float), but when I tried that on (complex single-float), it would return without error but with the result computed as if the imaginary part was 0. I thought maybe there was a padding issue, so I packed the two single-floats into a single 64-bit word that I then interpreted as a double-float. When passing this single argument to a function that expects a (complex single-float), the correct answer is computed.

For example, the BLAS function caxpy works on single-floats:

(FOREIGN-FUNCALL "gsl_blas_caxpy" 
		 :DOUBLE (PACK-COMPLEX-AS-DOUBLE ALPHA) 
		 :POINTER (MPOINTER X) :POINTER (MPOINTER Y)
		 :INT)

(FOREIGN-FUNCALL "gsl_blas_zaxpy" 
		 :DOUBLE (REALPART ALPHA) 
		 :DOUBLE (IMAGPART ALPHA) 
		 :POINTER (MPOINTER X) :POINTER (MPOINTER Y) 
		 :INT)

This has now been incorporated into GSLL. I assume it is not portable; I am on Debian amd64 but it works on x86 (32 bit) as well. It seems to be portable across CL implementations, at least it works in SBCL and CCL. So I have built in a test that tries a simple function; if it gets the wrong answer than all the functions that want a complex scalar will signal an error.

03:00 pm - GSLL new version

The new version GNU Scientific Library for Lisp (GSLL) is now available. This library has many mathematical functions used in science and engineering applications. It works on several Lisp implementations.

The most significant changes from the previous version include:

• All objects are memory managed (garbage collected), which means they can have indefinite extent like any other Lisp object. There is no need to put the object creation in a letm; in fact letm doesn't exist anymore. Use let, or defparameter, etc.
  
• All array (vector and matrix) element types that are supported by the platform, CL implementation, CFFI, and GSL, are supported by GSLL.
  
• On SBCL both Lisp and C use the same representation for array contents; they are not copied between the two sides. Thanks to Tamas Papp's foreign-friendly arrays for the inspiration.
  

There have been some function name and argument changes, so users of the previous version will need to update their code. Please use the git repository and discontinue using the old svn repository.

Feedback welcome, here or on the mailing list.

04:48 pm - Clozure CL on Debian

A while back, I mentioned trying OpenMCL for Debian. OpenMCL is now Clozure CL (CCL), and version 1.2 for MacOSX and Linux (PPC and amd64) was recently announced. I gave it another try and the experience is much better, though it needs Debianization.

The Debian Common Lisp Team has produced a Common Lisp in Debian manual, so I followed that to Debianize CCL. As root,

> lx86cl64
Welcome to Clozure Common Lisp Version 1.2-r10552  (LinuxX8664)!
? (load "/usr/share/common-lisp/source/common-lisp-controller/common-lisp-controller.lisp")
? (common-lisp-controller:compile-common-lisp-controller-v5 "ccl")
(#P"/var/cache/common-lisp-controller/0/ccl/common-lisp-controller/common-lisp-controller.lx64fsl"
 #P"/var/cache/common-lisp-controller/0/ccl/cl-asdf/asdf.lx64fsl"
 #P"/var/cache/common-lisp-controller/0/ccl/cl-asdf/wild-modules.lx64fsl"
 #P"/var/cache/common-lisp-controller/0/ccl/common-lisp-controller/post-sysdef-install.lx64fsl")

? (load "/var/cache/common-lisp-controller/0/ccl/common-lisp-controller/common-lisp-controller.lx64fsl")
? (load "/var/cache/common-lisp-controller/0/ccl/cl-asdf/asdf.lx64fsl")
? (load "/var/cache/common-lisp-controller/0/ccl/cl-asdf/wild-modules.lx64fsl")
? (load "/var/cache/common-lisp-controller/0/ccl/common-lisp-controller/post-sysdef-install.lx64fsl")
? (common-lisp-controller:init-common-lisp-controller-v5 "ccl")
? (ccl:save-application "ccl-deb")

I tried to compile and load various debian packages and other systems, with reasonably good success;

(clc:clc-require :asdf-system-connections)
(clc:clc-require :cl-base64)
(clc:clc-require :cl-ppcre)
(clc:clc-require :md5)
(clc:clc-require :port)
(clc:clc-require :iterate)
(clc:clc-require :cl-utilities)
(clc:clc-require :uffi)
(clc:clc-require :clsql-uffi)
(clc:clc-require :clsql-sqlite3)
(clc:clc-require :cffi)
(clc:clc-require :drakma)
(clc:clc-require :cl-html-parse)

(clc:clc-require :memoization)
 Error: The function CCL:ARGLIST is predefined in OpenMCL.
(use-package :iterate :cl-user)
> Error: Using #<Package "ITERATE"> in #<Package "COMMON-LISP-USER">
>        would cause name conflicts with symbols inherited by that package: 
>        ITERATE:TERMINATE  TERMINATE
 (clc:clc-require :trivial-http)
> Error: Module TRIVIAL-HTTP was not provided by any function on *MODULE-PROVIDER-FUNCTIONS*.

09:35 pm - SBCL in Debian testing (lenny): good news

Thanks to lots of work and attention from various people:

http://release.debian.org/migration/testing.pl?package=sbcl
Checking sbcl

* trying to update sbcl from 1:0.9.16.0-1 to 1:1.0.16.0-2 (candidate is 12 days old)
* sbcl is going in today

07:33 pm - SBCL in Debian testing (lenny)

Recently I noticed that SBCL in Debian testing is quite old; it is currently at 0.9.16 and evidently blocked there:

=== sbcl:
= Missing build(s) on alpha,sparc
 This might need manual action from your side.
 See http://buildd.debian.org/pkg.cgi?pkg=sbcl
= No migration to testing for 585 days.
 See <http://release.debian.org/migration/testing.pl?package=sbcl>

08:50 am - Cairo for Lisp

I have been looking into using Cairo for various 2D drawing needs. A key advantage is the availability of a variety of output formats, including both screen and hardcopy. One potential application is traditional x-y scientific/engineering plots. Although there are many options for doing that, even using Lisp, I have a Lisp package I adapted many years ago for this purpose that now lacks a good back end for output. My thought is that Cairo would be good for that, as it is a standardized library, widely available, with many modern output formats. I also have other 2D output that I do from Lisp, which I now do with a thin postscript definition layer that I wrote. Switching to Cairo would give me some more flexibility.

I did some research and found that there were three sets of bindings for Cairo:

• cl-cairo


• cl-cairo2


• cffi-cairo

It looks like I'll be using cffi-cairo.

02:16 pm - Never a NaN or Inf

In very few cases is it acceptable to me that a floating-point calculation would produce NaN or Inf; it almost always means that there is an error somewhere in my program. Following the principle that it is best to detect an error as soon as possible, I would like to know immediately whenever either of these is created, by trapping the exception and signaling an error. Depending on the language and implementation, this doesn't always happen by default. Here are some different languages and compilers, and how to set to make this happen:

• Lisp: SBCL: (sb-int:set-floating-point-modes :traps '(:invalid :divide-by-zero :overflow)) and
CMUCL: (ext:set-floating-point-modes :traps '(:invalid :divide-by-zero :overflow)).
  Also GSLL provides #'set-floating-point-modes.
  
• Fortran: for gfortran, use the compiler flag -ffpe-trap=invalid,zero,overflow.

• C/C++: use
  

  #include <fenv.h>
  trapfpe () {
    feenableexcept(FE_INVALID|FE_DIVBYZERO|FE_OVERFLOW);
  }
  

  
  In C, this will need to be included: void trapfpe();
  In C++, this will need to be included (also works for C):
  

  #ifndef _TRAPFPE_H_
  #define _TRAPFPE_H_
  
  #include <sys/cdefs.h>
  
  __BEGIN_DECLS
  void trapfpe(void);
  __END_DECLS
  
  #endif /* _TRAPFPE_H_ */
  

  
  and then in either, include trapfpe(); at the beginning. This should be portable to C99 compilers/libc but I've used only gcc and glibc. Also, GSL has gsl_ieee_env_setup.
  

08:54 am - SBCL, libc6, and GCC

For those that don't follow the mailing list, SBCL on Debian is unusable in unstable. SBCL did not change, but Debian upgraded libc6 to 2.7-9. This causes SBCL to spin at 100% usage, or simply to hang, while starting up. The only relevant change to libc6 is that it is now compiled with GCC 4.3 instead of 4.2; the libc6 maintainer has confirmed that backing out this change and recompiling libc6 with GCC 4.2 allows SBCL to work. It is still unclear where the problem(s) is/are; but the experts have narrowed it down further. Follow the developments on the Debian bug report; even if you don't use Debian, presumably this combination of software is a problem. In the meantime, if you are a sid user, don't upgrade! Though apparently some people have experienced the problem as a random failure, my experience is that it's 100% reproduceible.

Software is complicated.

06:27 pm - GNU Scientific Library for Lisp

The GNU Scientific Library is a library of applied mathemetics commonly used in science and engineering. I have written GNU Scientific Library for Lisp (GSLL), a fairly complete interface to this library from Common Lisp. My intent is that the interface be as Lisp-natural as possible. Though this will be useful for those of us that do scientific programming in Lisp, even those who aren't Lisp programmers might use this library as a desk-calculator interface to GSL.

01:29 pm - Comparison of floating point numbers in Common Lisp

Floating point numbers are a computer representation of the real numbers. Unlike the real numbers, there are a finite number of them. So there is a smallest and largest floating point number, and all others have a predecessor and successor.

Because different compilers and platforms can reorder a calculation and optimize in a way that is approximated differently and so do not necessarily produce the same floating point number, it is difficult to compare two floating point numbers and conclude that they represent the same result. For the purposes of regression (or unit) testing, we would like to do exactly this. Bruce Dawson addressed this problem in "Comparing Floating Point Numbers." He makes the point that the best way to do this correctly is to interpret each floating point number as an integer. By taking advantage of the IEEE 754 standard for representation of floating point numbers, we can construct a function that maps the floating point numbers to the integers. The genius of the standard's inventor W. Kahan is that a mapping derived from the standard, call it i(x), satisfies three properties:

• If two floats a<b, then i(a)<i(b),
• if two floats are adjacent and a<b, then i(b)=i(a)+1,
• and finally i(0.0)=0.

Dawson provides some clever C constructs to read a floating point number as an integer, and instruction on how to prevent the compiler from complaining about your trickery in doing so. Common Lisp instead provides us with functions with which we can properly construct our own integers. As a side benefit, we don't care what the actual representation of the floating point number is; we will build our own IEEE-like representation. We don't exactly want the full IEEE754 word though; we leave off the most significant bit, which is a sign bit, and instead make the sign of the integer agree with the sign of the float.

What we end up with is an enumeration of the floats. That is, for every single precision float, there is one integer in the range [-2139095039, 2139095039], and vice versa, with the exception that both positive and negative zero (allowed by the standard) map to zero. Likewise, there is a one-to-one mapping of the double precision floats to [-9218868437227405311,9218868437227405311]. Floats with special values (positive and negative infinity, and NaN) that are required by the standard do not have integer values.

I have written the following functions in Common Lisp:

• float-as-integer which is the function i(x);
• integer-as-float which is the inverse function (this isn't necessary but can be useful) and also returns the rational form of the float;
• decode-IEEE754 (used by other functions) that returns five values: significand, exponent, sign, bits in significand, bits in exponent, all as integers;
• format-IEEE754-bits which prints out the binary form of the IEEE word, separated into the three parts (this isn't necessary but is nice for comparing with bit expansions shown in references like the Wikipedia page).

Here are some interesting floats:

(float-as-integer most-negative-single-float)
-2139095039
(float-as-integer least-negative-single-float)
-1
(float-as-integer -0.0f0)
0
(float-as-integer 0.0f0)
0
(float-as-integer least-positive-single-float)
1
(float-as-integer (- 1.0f0 single-float-negative-epsilon))
1065353215
(float-as-integer 1.0f0)
1065353216
(float-as-integer (+ 1.0f0 single-float-epsilon))
1065353217
(float-as-integer most-positive-single-float)
2139095039

A regression test would record not the floating point number, but the integer produced by float-as-integer. Since integers can be unambiguously formatted to and read from a text file in a unique way, a subsequent recomputation would provide a clear indication of how close the floats are. Of course, we must decide how much error we're going to allow, because a correct calculation may produce slightly different integers. As an added bonus, these functions can be used to identify (in languages other than Lisp) when a positive single float has been interpreted as a double float.

06:00 pm - Multiprocessing lisp evaluations

I sometimes need to evaluate the same form with different parameters repeatedly. Such as

(job 1)
(job 2)
(job 3)
...

(defparameter *job-lock* (sb-thread:make-mutex :name "job lock"))
(defparameter *results* (list nil))
(defvar *end-of-jobs* (make-symbol "EOJ"))
(defvar *jobs* nil)

(defun worker (job)
  (let ((my-job nil)
	(more-jobs t))
    (loop while more-jobs
       do
       (sb-thread:with-mutex (*job-lock*)
	 (setf more-jobs (or *jobs*))
	 (setf my-job (when more-jobs (pop *jobs*))))
       (when my-job			; there is a job to be done
	 (if (eq my-job *end-of-jobs*)
	     (setf more-jobs nil)
	     (let ((my-results (apply job my-job))) ; call the job outside the mutex
	       (sb-thread:with-mutex (*job-lock*)   ; save results
		 (push my-results *results*))))))))

(defun run-tasks (job dataset number-of-workers)
  "The job is a function that takes one non-null argument.
   The dataset is a list of arglist sets for the job.
   The number-of-workers is the number of workers desired, 
   presumably the number of processors available."
  (setf *jobs* (make-list number-of-workers :initial-element *end-of-jobs*)
	*results* (list nil))
  (dolist (ds dataset) (push ds *jobs*))
  (let ((threads (list nil)))
    (loop repeat number-of-workers
       do (push (sb-thread:make-thread (lambda () (worker job))) threads))
    (dolist (thread (butlast threads)) (sb-thread:join-thread thread))
    (butlast *results*)))

(defun job (x) (list x (+ (loop for i from 1 to 2000000 sum (let ((p (* i x))) (1- (expt p (/ p))))))))
(run-tasks #'job '((1) (2) (3) (4) (5) (6) (7) (8) (9) (10)) 2)
((1 105.73587) (2 58.11023) (3 41.245914) (4 31.676012) (5 26.053244)
 (6 22.305136) (7 19.607342) (8 17.642727) (9 15.433696) (10 14.060191))

06:54 pm - Profiling in SLIME

This is about the easiest profiling I've seen in any language. In
fact, I think it's the only time I been able to make significant
improvements based on the report.

M-x slime-toggle-profile-fdefinition

M-x slime-profile-reset

M-x slime-profile-report

03:45 pm - Extension functions in sqlite3, again

I mentioned before that I cleaned up and posted mathematical and string extension functions for SQLite3. Some things have changed in the interim, which necessitated going through a few revisions. So now the new version is available which makes compilation and use easier. For one thing, the SQLite source is no longer required to compile it. Second, it uses the standard sqlite3_load_extension interface, which should make it easier to use. The interface through CL using CLSQL is now

(in-package :sqlite3)

;;; Add mathematical and string functions to SQL queries using
;;; libsqlitefunctions from
;;; http://www.sqlite.org/contrib/download/extension-functions.c?get=22

(def-sqlite3-function
    "sqlite3_enable_load_extension"
    ((db sqlite3-db) (onoff :int))
 :returning :int)

(def-sqlite3-function
    "sqlite3_load_extension"
    ((db sqlite3-db)
     (filename :cstring)
     (entrypoint :int)
     (errmsg :int))
  :returning :int)

(eval-when  (:compile-toplevel :load-toplevel :execute)
  (export 'enable-sqlite3-extension-functions))

(defun enable-sqlite3-extension-functions (database)
  "Set up the SQLite3 mathematical extension functions.  This
   must be called every time the database is connected
   before any extension function is used."
  (let ((db-ptr (clsql-sqlite3::sqlite3-db database)))
    (sqlite3-enable-load-extension db-ptr 1)
    (unless (zerop
	     (sqlite3-load-extension db-ptr "libsqlitefunctions.so" 0 0))
      (error "Can't load libsqlitefunctions.so."))))

11:57 am - cl-mcclim out of date in Debian unstable

People are talking about the hot new mcclim (0.9.5). But unstable has 0.9.2; experimental only has 0.9.4 and it has not migrated to unstable in the six months it's been there because, supposedly, the etch freeze http://packages.qa.debian.org/c/cl-mcclim.html. But etch has long since been released; why the continued logjam?

01:39 pm - DRAKMA

Yesterday's blog mentioned trouble I'm having with cl-curl, my own package. A comment pointed me to DRAKMA, yet another Edi Weitz bequest to the lisp world. It is a native common lisp http client that handles the particular web site I am scraping data from, complete with password and cookies. It was very simple to rewrite my access functions to use DRAKMA, and it is able to retrieve all the data I wanted, no memory fault like cl-curl. So I will update the cl-curl project page to recommend DRAKMA.

It was a little unclear how to install DRAKMA; here is my condensed summary, on Debian:

1. sudo aptitude install cl-chunga cl-puri cl-flexi-streams
   
2. wget http://common-lisp.net/project/usocket/releases/usocket-0.3.2.tar.gz
   
3. wget http://common-lisp.net/project/cl-plus-ssl/download/cl+ssl-2007-03-10.tar.gz
   
4. wget http://weitz.de/files/drakma.tar.gz
   
5. tar zxvf usocket-0.3.2.tar.gz
   
6. tar zxvf cl+ssl-2007-03-10.tar.gz
   
7. tar zxvf drakma.tar.gz
   
8. clc-register-user-package usocket-0.3.2/usocket.asd
   
9. clc-register-user-package cl+ssl-2007-03-10/cl+ssl.asd
   
10. clc-register-user-package drakma-0.7.0/drakma.asd
    
11. ,l drakma (in slime)
    

10:24 am - Parsing HTML and memory fault from cl-curl

I've had need to parse HTML in lisp from time to time. The latest reason is some very specific and uncomplicated HTML that I scrape for some satellite data off a published database. A search online turns up XMLS as a likely candidate for this task. I have used it successfully in the past, but recently I find it won't parse its own example and its own supplied HTML documentation, to say nothing of the real HTML I want it to parse. It either returns NIL (meaning an error in parsing the correct HTML), or only the first line of HTML. There is a thread on comp.lang.lisp about how to parse HTML, and many people recommend cl-html-parse. I was dissuaded at first because of the wiki comments implying that it had been superseded by pxmlutils whose web page in turn implies that it had been superseded by... XMLS! But cl-html-parse works just fine on the web pages I need to scrape.

So, success. But then, I am grabbing the web page with cl-curl which works most of the time, but for a particular query gives "memory fault," I think because there is a lot of data. And the author/maintainer of cl-curl is... me! D'Oh. It would be nice to have cl-curl using CFFI instead of UFFI, maybe based on the pedagogical development given in the CFFI tutorial. My hope is at least it would solve this problem. Any interest/volunteers/motivators?

09:59 am - SQLite math and string functions simplified

Thanks to Mikey C, there are now a number of math and string functions available in SQLite SQL queries. I blogged this before, but now I have simplified and cleaned up his code, and made it available.

The functions are: Math: acos, asin, atan, atn2, atan2, acosh, asinh, atanh, difference, degrees, radians, cos, sin, tan, cot, cosh, sinh, tanh, coth, exp, log, log10, power, sign, sqrt, square, ceil, floor, pi. String: replicate, charindex, leftstr, rightstr, ltrim, rtrim, trim, replace, reverse, proper, padl, padr, padc, strfilter.

For CLSQL users, the Lisp portion of the instructions in my previous post remain the same.