ANNOUNCE: normaldistribution-1.1 – Minimum fuss normally distributed random values.

I’m pleased to announce the immediate availability of the normaldistribution library on hackage.

This purpose of this library is to have a simple API and no dependencies beyond Haskell 98 in order to let you produce normally distributed random values with a minimum of fuss. This library does not attempt to be blazingly fast nor to pass stringent tests of randomness. It attempts to be very easy to install and use while being “good enough” for many applications (simulations, games, etc.). The API builds upon and is largely analogous to that of the Haskell 98 Random module (more recently System.Random).

Pure:

> (sample,g) = normal  myRandomGen  -- using a Random.RandomGen
> samples    = normals myRandomGen  -- infinite list
> samples2   = mkNormals 10831452   -- infinite list using a seed

In the IO monad:

> sample    <- normalIO
> samples   <- normalsIO  -- infinite list

With custom mean and standard deviation:

> (sample,g) = normal'    (mean,sigma) myRandomGen
> samples    = normals'   (mean,sigma) myRandomGen
> samples2   = mkNormals' (mean,sigma) 10831452

> sample    <- normalIO'  (mean,sigma)
> samples   <- normalsIO' (mean,sigma)

Internally the library uses the Box-Muller method to generate normally distributed values from uniformly distributed random values. If more than one sample is needed taking samples off an infinite list (created by e.g. normals) will be roughly twice as efficient as repetedly generating individual samples with e.g. normal.

Update 2011–04–09: Changed to reflect version 1.1.

Update 2011–04–22: Version 1.1.0.1 released which builds with haskell98-1.1+ (GHC-7.0.1+). No code changes.

Flattrable Reactive Programming

I have occasionally seen these “flattr this!” buttons, notably on the blog of Luke Palmer, and occasionally felt it was exactly what I wanted to do. However, inertia always got the better of me – no doubt I was busy being distracted by some lolcat screaming for adoration.

Lately I noticed a flattr button on the blog of Heinrich Apfelmus too, and again felt a desire to give. Heinrich’s blog was the tipping point (pun intended) that finally got me to actually visit flattr.com, and this great quote from Alexandros closed the sale:

This month, thanks to Flattr, I’we shared my world to the world. I earned 0.03 cents - the best 0.03 cents ever!

How can that not make me want to give? (Is Alexandros being sarcastic? He looks happy enough in his profile pic…)

My first flattrs go to Luke and Heinrich; 0.03 cents are coming your way! ;)

Alexandros gets flattred too. And now that I’ve signed up I’ll be on the lookout for more flattr buttons (so don’t feel bad if you weren’t mentioned here!).

As an aside, the common denominator that made me want to flattr Luke and Heinrich is their efforts to explore Functional Reactive Programming. It seems to me that much remains to be discovered in the FRP design space and I hope my flattrs encourage these brave explorers to push the boundaries further… Oh, look! Conal Elliott has a flattr button too!

Dimensional example: leaky container pressure

Herein I demonstrate a calculation performed today in my work. The details are slightly different to protect the innocent but the problem is the same. The purpose of all this is to demonstrate use of the dimensional library. What this blog post unfortunately fails to demonstrate is how the compiler helps catch errors in units and formulae, effectively double-checking any derivations for me – this is the real value of dimensional! In case you have never heard of dimensional before here is the blurb from the project site:

Dimensional is a library providing data types for performing arithmetic with physical quantities and units. Information about the physical dimensions of the quantities/units is embedded in their types and the validity of operations is verified by the type checker at compile time. The boxing and unboxing of numerical values as quantities is done by multiplication and division with units. The library is designed to, as far as is practical, enforce/encourage best practices of unit usage.

This is a literate Haskell program. You can copy and paste this blog post into a .lhs file and it should compile and run (provided you have installed dimensional). First the formalities:

> import Numeric.Units.Dimensional.Prelude
> import Numeric.Units.Dimensional.NonSI
> import qualified Prelude

(In the below I have used parentheses rather than juxtaposition to denote function application for functions of time, e.g. p(t). Does this help or hurt clarity given the subject matter?)

Consider a leaky container of volume v filled with gas to pressure p₀ and launched into space at time t₀ (insert favorite monad tutorial joke here). There is a requirement that at some future time t₁ said container must retain no less than a given amount of gas. The temperature T of the container is held constant by a thermal control system so per the ideal gas law (a good enough approximation) we can characterize the requirement to remnant gas by its pressure p_min. The problem at hand was to determine the maximum allowable leak rate of the container at t₀ to ensure that the pressure at time t₁ is no less than p_min.

Here are the inputs (don’t blame me for the choice of units, this is what I was given):

> v     =      2 *~ liter
> t_0   =      0 *~ hour
> p_0   =    750 *~ mmHg  -- roughly atmospheric pressure
> t_1   = 133650 *~ hour  -- about 15.25 years
> p_min =      5 *~ mmHg
> temp  = fromDegreeCelsiusAbsolute 22

For slow leakage it can be assumed that the leakage rate is proportional to the pressure of the gas inside the container (more generally, the pressure difference between the gas in the container and the surrounding medium). Thus the rate of change of the pressure is described by the first order linear differential equation

> dpdt(t) = negate k * p(t)

with solution

> p(t) = c * exp (negate k * t)

The integration constant c is determined from the initial conditions at t₀:

> c = p_0  -- follows from t_0 = 0 s.

We chose the worst case conditions at t₁ to determine k:

> k = log (c / p_min) / t_1

Now, per the ideal gas law the amount of substance n in the container is described by

> n(t) = p(t) * v / (r * temp)

where

> r = 8.314472 *~ (joule / (mole * kelvin))

is the ideal gas constant. Differentiating n(t) give the rate of change of substance in the container:

> dndt(t) = dpdt(t) * v / (r * temp)

and the maximum allowed rate of change (leakage rate) at time t₀ is produced by

> main = print (dndt t_0)  -- -8.486687441280605e-10 s^-1 mol

That’s it!

ANNOUNCE: numtype 1.0 — Type-level (low cardinality) integers

Since its inception my dimensional library has been built around a unary type-level representation of integers (NumTypes) defined in the Numeric.NumType module. This module has proven itself useful outside the context of dimensional and after dragging my feet for a long time I’ve finally gotten around packaging it up in its own library: numtype.

The Numeric.NumType module is completely self-contained (only imports Prelude) and is heavily commented in a narrative manner inspired by Oleg Kiselyov’s expositions. I believe it provides a good case study for type-level programming with multi-parameter type classes (MPTCs) and functional dependencies.

Addition, subtraction, division, and multiplication of NumTypes is supported. NumTypes have no value-level representation but can be converted to any Num instance with the toNum function.

The numtype library has two significant short-comings:

• Minimal haddocks — as with my dimensional library the literate Haskell source code is the documentation. The flip-side is that the code is very well-commented.

• Due to the unary implementation the practical size of the NumTypes is severely limited, making them unsuitable for large-cardinality applications. If you will be working with integers beyond (–20, 20) this package probably isn’t for you.

(If the second bullet is a show-stopper Edward Kmett’s type-int library may be a better choice. Peter Gavin’s tfp library also provides type-level integers but uses type families instead of MPTCs and fundeps. I cannot vouch for either of these libraries as I haven’t used them.)

Numtype version 1.0 can be downloaded from Hackage or the dimensional project page. I’ve also updated dimensional to version 0.8 with the Numeric.NumType module removed and Julian year and century units added. Enjoy!

Benchmarking Amazon EC2 with GHC

My personal computers are pretty old and/or slow. I have an old PowerBook G4 and a newish EEE PC. The PowerBook was top of its class when I got it with most options maxed out. Alas, that was five years ago. The EEE PC is by definition not a top performer, nor does it try to be. I find that when the two machines perform similarly in day-to-day tasks, at least when the EEE PC is in “Super performance” mode.

Truth is that for most of the mundane stuff I do these two machines perform acceptably. I’m obviously not going to be watching any 1080p movies on them, or enjoying the latest games (from a productivity standpoint I’m not sure these are such bad things), but most everything else works fine.

Where I feel the performance does hurt me is when compiling with GHC (or any compiler for that matter, it’s just that GHC is the one I use most). Often I spend way to much of my precious little private developer time waiting for the compiler to finish. This is in stark contrast to the situation at work where I run GHC on a pretty zippy Dell PowerEdge Blade Server.

In the near future I expect to be spending more time developing at home and want to be able to do so more efficiently, preferably on par with the situation at work. I could obviously buy shiny new hardware but being a miser (in case you couldn’t already tell based on my hardware) I’m looking for alternatives that would allow me to avoid or defer a hefty up-front investment. One such alternative I’m considering is to rent compute capacity in the Amazon Elastic Compute Cloud (EC2).

EC2 compute capacity is sold in the form of instances at an hourly rate ranging from $0.10 to $0.80 depending on capacity/performance plus some small-change for bandwidth and persistent storage. While there are a couple of hurdles to overcome in order to leverage EC2 as a development workstation the first thing I want to do is to make sure it is a goal worthy of pursuing in the first place, i.e. will I get the desired performance gains at a reasonable price?

To at least begin to answer this question I’ve done some informal benchmarking of the aforementioned systems, excluding the pricier EC2 instances. All the regular benchmarking caveats apply and to reinforce the unscientificity of it all I’m not going to bother providing complete specs for the systems. Here are the fundamentals:

• Apple PowerBook G4: 1.5 GHz PowerPC G4 processor, 2 GB RAM, 5400 rpm HD.
• Asus EEE PC 900HA: 1.6 GHz Intel Atom processor, 1 GB RAM, 4200 rpm HD.
• Dell PowerEdge 1855 Blade Server: Two single-core 3.2 GHz Xeon processors, 2 GB RAM.
• Amazon EC2 Small Instance: 1 EC2 Compute Unit (1 virtual core), 1.7 GB RAM, $0.10 per hour.
• Amazon EC2 High-CPU Medium Instance: 5 EC2 Compute Units (2 virtual cores), 1.7 GB RAM, $0.20 per hour.

According to Amazon “one EC2 Compute Unit provides the equivalent CPU capacity of a 1.0–1.2 GHz 2007 Opteron or 2007 Xeon processor.”

I figure the results are probably more interesting that the details of the tests so here they are, the systems are ordered by increasing performance which happened to be consistent across the tests:

Benchmark results, times in seconds, shorter times are better.

	astro-tables build	highlighting-kate build	fad test suite
PowerBook G4	292	(848)	28
EEE PC	291	643	18
EC2 Small	171	519	15
Dell PowerEdge	75	260	6
EC2 Medium	55	172	4

The EEE PC was in “Super Performance” mode during the tests and the PowerBook was at its highest CPU speed. All times are the “real” time as measured by the Unix time command and lower numbers are naturally better.

As can be seen an Amazon EC2 Small instance is only marginally faster than the EEE PC. An EC2 High-CPU Medium instance on the other hand is significantly faster than the zippy Dell PowerEdge Blade server. Is either one a good deal? Good question, I think a case could be made either way depending on your priorities but I’m not going to tackle that today.

If you care about the details of the tests read on, if not please move on to your next blog of choice!

astro-tables build

This package currently consists of a single automatically generated 4000-line monster of a module¹. The code is pretty straight-forward but the module takes ages to compile, almost certainly due to me giving the type checker an unnecessarily hard time. I have a trivial rewrite on my todo-list which I expect will shorten the compilation time dramatically but the current form comes in kind of handy for the purposes of this benchmark. The Git repo is git://github.com/bjornbm/astro-tables.git and the commit used in the benchmarking was e63b8978833878526870b2101697197ff64af593. I made sure the dependencies were already installed and ran time cabal install.

highlighting-kate build

From recent memory I knew that John MacFarlane’s highlighting-kate package has a hefty number of modules (the majority of which are also automatically generated) that take a fair amount of time to compile. I downloaded version 0.2.4 from hackage, made sure all dependencies were already installed, and ran time cabal install --flags=executable.

I ran into one snag with this test: the build wouldn’t complete with GHC 6.10.3 on the PowerBook G4 due to some problem with pcre-light (which tends to give me headaches on pretty much every platform). This particular headache² I was unable to resolve and had to run the test using GHC 6.10.1 on the PowerBook G4.

fad test suite

Finally I did a runtime performance (as opposed to compilation) benchmark: running the test suite of the fad library. The Git repo is git://github.com/bjornbm/fad.git and the commit used was cd2965a6741291570930e4bf6e9f8f9ab64ccadd. I ran ghc --make Test and then time ./Test.

---

1. An implementation of the 678 lunisolar terms and 687 planetary terms of the IAU 2000A Precession-Nutation Model. ↩

2. gcc: Internal error: Virtual timer expired (program cc1) ↩

May 2009 HCAR Submissions

Last weekend was the submission deadline for contributions to the May 2009 edition of the Haskell Communities and Activities Report (HCAR). I absolutely love the HCAR; it’s a terrific and comprehensive source of information on all the amazing stuff people are doing with Haskell. I’m very grateful to the editor Janis Voightländer (and Andres Löh before him) for compiling the report, in my opinion he does a great service to the Haskell community. I can’t wait to get my hands on the upcoming edition!

Below are my submissions to the HCAR (free blogging material, yay!) with some additional notes.

dimensional: Statically checked physical dimensions

Report by: Björn Buckwalter
Status: active, mostly stable

Dimensional is a library providing data types for performing arithmetics with physical quantities and units. Information about the physical dimensions of the quantities/units is embedded in their types, and the validity of operations is verified by the type checker at compile time. The boxing and unboxing of numerical values as quantities is done by multiplication and division with units. The library is designed to, as far as is practical, enforce/encourage best practices of unit usage.

The core of dimensional is stable with additional units being added on an as-needed basis. In addition to the SI system of units, dimensional has experimental support for user-defined dimensions and a proof-of-concept implementation of the CGS system of units. I am also experimenting with forward automatic differentiation and rudimentary linear algebra.

The current release is compatible with GHC 6.6.x and above and can be downloaded from Hackage or the project web site. The primary documentation is the literate Haskell source code, but the wiki on the project web site has a few usage examples to help with getting started.

Further Reading

http://dimensional.googlecode.com

Dimensional was largely the project that enticed me to learn Haskell, and is the one that I am most proud of to date. It relies on some pretty tricky type-level hackery, stuff I would never have been able to pull off without the help of many papers ranging from Oleg Kiselyov’s amazing type-hacks to Hudak et al’s A History of Haskell.

Its primary shortcoming at present is the lack of Haddock documentation. Dimensional is written in a literate style and I think anyone reading the source will find the documentation pretty comprehensive and easy to follow (barring the fact that type hacking is pretty tricky business). Alas, with the last few years’ infrastructure improvements and increased (de facto) standardization in the community not having proper haddocks is a wart.

fad: Forward Automatic Differentiation

Report by: Björn Buckwalter
Participants: Barak A. Pearlmutter, Jeffrey Mark Siskind
Status: active

Fad is an attempt to make as comprehensive and usable a forward automatic differentiation (AD) library as is possible in Haskell. Fad (a) attempts to be correct, by making it difficult to accidentally get a numerically incorrect derivative; (b) item provides not only first-derivatives, but also a lazy tower of higher-order derivatives; (c) allows nested use of derivative operators while using the type system to reject incorrect nesting (perturbation confusion); (d) attempts to be complete, in the sense of allowing calculation of derivatives of functions defined using a large variety of Haskell constructs; and (e) tries to be efficient, in the sense of both the defining properties of forward automatic differentiation and in keeping the constant factor overhead as low as possible.

Version 1.0 of fad was uploaded to Hackage on April 3. Recent changes can be found via git clone git://github.com/bjornbm/fad.git

Further Reading

http://github.com/bjornbm/fad
http://flygdynamikern.blogspot.com/2009/04/announce-fad–10-forward-automatic.html

The majority of work on fad (including most of the above report) is actually being done by Barak Pearlmutter, principial of the Hamilton Institute’s Brain and Computation Lab. As you can understand his brain is in a different league than mine when it comes to anything related to computation. My contribution at this stage is mostly reviewing his patches (often struggling to understand them) and providing guidance and support on Haskell infrastructure and conventions.

leapseconds-announced

Report by: Björn Buckwalter
Status: stable, maintained

The leapseconds-announced library provides an easy to use static LeapSecondTable with the leap seconds announced at library release time. It is intended as a quick-and-dirty leap second solution for one-off analyses concerned only with the past and present (i.e. up until the next as of yet unannounced leap second), or for applications which can afford to be recompiled against an updated library as often as every six months.

Version 2009 of leapseconds-announced contains all leap seconds up to 2009–01–01. A new version will be uploaded if/when the IERS announces a new leap second.

Further Reading

http://hackage.haskell.org/cgi-bin/hackage-scripts/package/leapseconds-announced
http://github.com/bjornbm/leapseconds-announced

Not much to add about leapseconds-announced other than that it fills a need I have. I elaborated a little on what it is (and what it doesn’t try to be) in the announcement thread.

That’s it! I hope to have another project or two to report on in the November edition…

ANNOUNCE: fad 1.0 — Forward Automatic Differentiation for Haskell

I’m pleased to announce the initial release of the Haskell fad library, developed by Barak A. Pearlmutter and Jeffrey Mark Siskind. Fad provides Forward Automatic Differentiation (AD) for functions polymorphic over instances of Num. There have been many Haskell implementations of forward AD, with varying levels of completeness, published in papers and blog posts, but alarmingly few of these have made it into hackage — to date Conal Elliott’s vector-spaces package is the only one I am aware of.

Fad is an attempt to make as comprehensive and usable a forward AD package as is possible in Haskell. However, correctness is given priority over ease of use, and this is in my opinion the defining quality of fad. Specifically, Fad leverages Haskell’s expressive type system to tackle the problem of perturbation confusion, brought to light in Pearlmutter and Siskind’s 2005 paper Perturbation Confusion and Referential Transparency. Fad prevents perturbation confusion by employing type-level “branding” as proposed by myself in a 2007 post to haskell-cafe. To the best of our knowledge all other forward AD implementations in Haskell are susceptible to perturbation confusion.

As this library has been in the works for quite some time it is worth noting that it hasn’t benefited from Conal’s ground-breaking work in the area. Once we wrap our heads around his beautiful constructs perhaps we’ll be able to borrow some tricks from him.

As mentioned already, fad was developed primarily by Barak A. Pearlmutter and Jeffrey Mark Siskind. My own contribution has been providing Haskell infrastructure support and wrapping up loose ends in order to get the library into a releasable state. Many thanks to Barak and Jeffrey for permitting me to release fad under the BSD license.

Fad resides on GitHub and hackage and is only a cabal install fad away! What follows is Fad’s README, refer to the haddocks for detailed documentation.

   Copyright  : 2008-2009, Barak A. Pearlmutter and Jeffrey Mark Siskind
   License    : BSD3

   Maintainer : bjorn.buckwalter@gmail.com
   Stability  : experimental
   Portability: GHC only?

Forward Automatic Differentiation via overloading to perform
nonstandard interpretation that replaces original numeric type with
corresponding generalized dual number type.

Each invocation of the differentiation function introduces a
distinct perturbation, which requires a distinct dual number type.
In order to prevent these from being confused, tagging, called
branding in the Haskell community, is used.  This seems to prevent
perturbation confusion, although it would be nice to have an actual
proof of this.  The technique does require adding invocations of
lift at appropriate places when nesting is present.

For more information on perturbation confusion and the solution
employed in this library see:
<http://www.bcl.hamilton.ie/~barak/papers/ifl2005.pdf>
<http://thread.gmane.org/gmane.comp.lang.haskell.cafe/22308/>


Installation
============
To install:
    cabal install

Or:
    runhaskell Setup.lhs configure
    runhaskell Setup.lhs build
    runhaskell Setup.lhs install


Examples
========
Define an example function 'f':

> import Numeric.FAD
> f x = 6 - 5 * x + x ^ 2  -- Our example function

Basic usage of the differentiation operator:

> y   = f 2              -- f(2) = 0
> y'  = diff f 2         -- First derivative f'(2) = -1
> y'' = diff (diff f) 2  -- Second derivative f''(2) = 2

List of derivatives:

> ys = take 3 $ diffs f 2  -- [0, -1, 2]

Example optimization method; find a zero using Newton's method:

> y_newton1 = zeroNewton f 0   -- converges to first zero at 2.0.
> y_newton2 = zeroNewton f 10  -- converges to second zero at 3.0.


Credits
=======
Authors: Copyright 2008,
Barak A. Pearlmutter <barak@cs.nuim.ie> &
Jeffrey Mark Siskind <qobi@purdue.edu>

Work started as stripped-down version of higher-order tower code
published by Jerzy Karczmarczuk <jerzy.karczmarczuk@info.unicaen.fr>
which used a non-standard standard prelude.

Initial perturbation-confusing code is a modified version of
<http://cdsmith.wordpress.com/2007/11/29/some-playing-with-derivatives/>

Tag trick, called "branding" in the Haskell community, from
Bjorn Buckwalter <bjorn.buckwalter@gmail.com>
<http://thread.gmane.org/gmane.comp.lang.haskell.cafe/22308/>