Class 12: Foreign Functions and Array Operations

Review from previous class

Moving average smoothers

These are the simplest to define and hint at the problems that need to be handled, such as the choice of bandwidth and the adequacy of a single bandwidth.

Foreign functions

External functions can be very helpful in making critical parts of your program much, much faster. How much depends on the application, but factors of 10 are not uncommon.

Status of Projects

• Need to start investing more of your time.
• My page is still ahead of most. That's not too good.

Importing Foreign Functions into Lisp, Part 2

Better smoothing code

As noted last time, its generally better to think harder about the algorithm before investing the time to write C code. At least try the compiler. Here's a much improved version of our moving average smoother that uses the smooth at time t S_t to compute S_t+1.

		   (defun MOVING-AVERAGE (x width)
		     "Smooths assuming input x is equally spaced."
		     (unless (oddp width)
		             (format t "Width needs to be odd integer >= 3.~%")
		             (setf width (let ((rw (round width))  )
		                           (if (oddp rw) rw (1+ rw))))   )
		     (let* ((half  (/ (1- width) 2))             ; invariants done once
		            (n     (length x))
		            (count half)
		            (sum   (sum (select x (iseq half)))) ; initialize terms
		            (drop  (- (1+ half)))
		            (add   half)
		            (smth  ())  )
		       (dotimes (i n)
		                (when (<= 0 drop)
		                      (decf sum (select x drop))
		                      (decf count))
		                (when (<  add n)
		                      (incf sum (select x add))
		                      (incf count))
		                (incf drop)
		                (incf add)
		                ; (format t "~d ~d  count=~d~%" drop add count)
		                (push (/ sum count) smth))
		       (nreverse smth)))
        

The speed of the calculation does not depend on the width of the smoother with this approach. However, as we have made it faster, we have also made it more specialized. This version does not generalize to weighted averages so simply as the prior version. We have built the constant smoothing weights into the algorithm, replacing the call to mean by the decf's and incf's.

In this form, animation works rather nicely as long as the series to be smoothed is no more than 100-200 points.

Foreign function version

Here's the C version of the same program. Note the obvious similarity of the code. The syntax differs here and there (for example, ++add rather than incf add) but the heart of the algorithm remains the same.

			void MovAvg (double x[], long *n, long *w, double smth[])
			{
				long i, half, count, drop, add;
				double sum;
				
				count = half = (*w - 1)/2;
				sum = 0.0;
				for (i=0; i < half; ++i)
					sum += x[i];
				drop = - (1+half);
				add = half;
				for (i=0; i< *n; ++i)
				{	if (drop >= 0)
					{	sum -= x[drop];
						--count;
					}
					if (add < *n)
					{	sum += x[add];
						++count;
					}
					smth[i] = sum / count;
					++drop;
					++add;
				}
			}
		

In order to use this function, we first compile it into an "object file" movAvg.o using the Unix cc command:

            cc -c +O3 -Aa +Z  movAvg.c
        

The man pages describe these options which

1. suspend the usual building of an application (-c),
2. control optimization (+O3),
3. describe the compiler options, here to use ANSI C (-Aa), and
4. make the code position independent (+Z).

Before we can load this code into the Lisp environment, we also need to wrap it up as a shared library so that the Lisp routines are able to extract the function references. The next command coverts the output file movAvg.o into a shared library called movAvg.shlb.

            ld -o movAvg.shlb -b movAvg.o
		

Now we can use dyn-load to read the shared library, and then use call-cfun to make things happen.

Arrays and Matrices in LispStat

Like common Lisp, LispStat includes basic features for handling arrays. Unlike standard common Lisp, thankfully, LispStat has a collection of helpful matrix algorithms which are important in statistics problems.

Lisp script for today's class

class12.lsp

Next time