Monday, September 17, 2018

New publications

Peer reviewed


"PC proxy: A method for dynamical tracer reconstruction," Environmental Fluid Mechanics, doi:10.1007/s10652-018-9615-7.

"Accelerating kernel classifiers through borders mapping," Journal of Real-Time Image Processing, doi:10.1007/s11554-018-0769-9.

Pre-prints/under review


"Solving for multi-class using orthogonal coding matrices," arXiv: 1801.09055.

"Solving for multi-class: a survey and synthesis," arXiv: 1809.05929.


A full list of my publications can be found on Google scholar. Note that the peer-reviewed versions of the first two papers differ from the pre-print submissions. In particular that of the second paper is significantly different from the final, published version. A one year embargo prevents me from making the published versions public, however.

Tuesday, July 4, 2017

How to write a scientific paper

First you need to come up with a stack of references. It always looks good to reviewers if you've got a reference list as long as your arm. You need a small army of other scientists on whose broad shoulders your authority may rest. It looks bad if each and every idea, no matter how trivial, is not at least hinted at in some other work. Heaven forbid you should come up with something original.

Always pepper your writing with plenty of faux Latin phrases. Terms like "in situ" "et al. "i.e." "vis a vis" "a priori" "post hoc" etc. It makes you sound educated.

Whenever you come up with a new method, algorithm or computer program, Christen it with a cute, catchy or flashy acronym. ARTS (Atmospheric Radiative Transfer Simulator) and ARTIST (I forget what this one stands for but it's for sea ice) are two nice ones from my own field. It's also good to abbreviate every technical phrase you use more than twice. It makes you look like you've got a firm, insider grasp of topics that no one else understands.

Try to squeeze in as many equations as possible. Concepts that are easy to explain in words can always be made more obtuse by converting them into formulas, er, sorry, formulae. Besides, like a long reference list, it looks impressive. If your reviewers don't have to take at least 3 Tylenol IIIs (six regular) after puzzling unsuccessfully over your article, it just isn't ready for prime time.

Sunday, April 3, 2016

The best derailleur ever made


The Huret Allvit? Really? Now hear me out. I know it has a bad reputation, but I'm not convinced that reputation is deserved. While the execution is not perfect, the design is positively inspired and it has a number of other design elements which, to me at least, scream "QUALITY". The problem with these design elements, and I will get into this, is that they are not idiot-proof. Consider:

1. The parallelogram design moves the top pulley further out as you move to larger sprockets in a design that predates the patented Suntour slant parallelogram by several years (1958 vs. 1964). Moreover, this motion is nonlinear--the motion is faster as the pulley cage is moved inwards, similar to how there is a bigger difference in teeth as the sprockets get larger.


2. The jockey pulleys have no teeth and fully adjustable ball bearings. Since the pulleys transmit no torque, there is no need for teeth. Both features will reduce drivetrain friction.

3. Every single pivot except for the cable hanger (which is riveted on) has a fully adjustable bearing with a locknut. Compare this with modern derailleurs which are almost always riveted together. Once the pivots wear out, the unit is done: you have to throw it away. With the Allvit, if the shifting becomes less precise because the pivots are loose, you just pull the thing off, disassemble it, clean up each part, and put it back together with the pivots nice and tight. This is almost certainly the source of the bad reputation these derailleurs had: if just one of these pivots is out of adjustment then it won't work well, if at all.

It was not all good, of course. The cable clamp was particularly poor and in conjunction with the strong parallelogram return spring and the thin gauge shifter cables of the time produced a lot of broken cables.

Towards the end of its life, the Allvit was being used on a lot of cheap, department store "ten-speed" and utility bikes. In order for it to work well, it needs to be maintained, something that's not going to happen with these cheaper bikes. The Japanese derailleurs, by contrast, are more idiot proof. They always work, no matter who attached them to the bike or how badly. They also tend to be more aesthetically pleasing. Note how on this vintage Shimano 600 derailleur the adjustment screws are placed neatly beside each other alongside the top parallelogram plate:
Contrast the more industrial look of the Allvit with one adjustment screw on the fixed top plate and the other down below jutting out of the bottom parallelogram arm.

Sunday, July 26, 2015

Move to github

I guess I've succumbed to pressure and moved over to github. Sourceforge was down for over a week and there has been a lot of bad press about it lately. I've rolled all of my software projects of any importance into a single scientific mega-library which you can find here:

libmsci

Not sure if this is a good approach. Probably not. I was just getting sick of dependency issues especially during release time. None of the libraries are all that big, so I thought, why not just roll them into one?

Also, be sure to check out the new homepage:

Peteysoft.github.io

Monday, June 1, 2015

Test driven development

Although I don't consider myself a computer programmer no matter how much code I churn out, nonetheless I feel I ought to keep up with some of the latest trends in software development.  Don't get me wrong, but I also suspect that many of these trends are also a lot of hot air and that they go in and out of fashion.  Certainly object-oriented programming is over-rated.

One idea that's caught my attention is so-called test-driven development.  I know that I don't write enough tests for my computer programs.

Unlike some other current fashions such as functional programming, this is something I can put into practice right away.  Again, I have my doubts about it so I will describe my initial attempts here.  I guess this is a somewhat boring topic, delving once again into the minutiae of scientific programming--hopefully I will have the wherewithal to write about more far-reaching and interesting topics later on some point.

To start with, I've picked an object class that's easy to test.  I have a family of classes that select out the k-least elements in a list.  Because some methods, such as a tree or heap, don't require storing the whole list, the object works by adding each item element-by-element and then extracting the k-least.  The base class looks like this:

template
class kiselect_base {
  protected:
    long ncur;
    long k;
  public:
    kiselect_base();
    kiselect_base(long k1);
    virtual ~kiselect_base();
    //add an element, set associated index to element count; 
    //return current size of data structure:
    virtual long add(type val)=0;
    //add an element and set the index, returns current size of data structure:
    virtual long add(type val, long ind)=0;
    //marshall out the k-least and associated indices:
    virtual void get(type * kleast, long *ind)=0;
    //test the selection algo:
    int test(long n); //number of elements
};

I wrote a family of these things because I wanted to test out which version was fastest.  It turns out that it makes fat little difference and even a naive method based on insertion or selection with supposed O(nk) performance is almost as good as a more sophisticated method based on quick-sort with supposed O(n) performance.  In addition to the k-least elements of the list, this version returns a set of numbers that index into the list in case there is auxiliary data that also needs to be selected.  This also makes slightly easier to test.

As you can see I've already added the test method.  Sticking it in the base class means that all of the children can be tested without any new code.  This is precisely in keeping with my thoughts about test routines: they should be as general as possible.  Forget having a small selection of specific test cases: this is a recipe for disaster.  Maybe it's not likely, but suppose your code just happens to pass all of them while still being incorrect?  Plus it's trivial to write code that passes all the test cases without being anywhere close to correct.

Rather, we need to be able to generate as many random test cases as we need.  Five test cases not enough?  How about 100?  How about one million?  Here is my first crack at the problem:

//trying to move towards more of a test-driven development:
template
int kiselect_base::test(long n) {
  int err=0;
  type list[n];
  type kleast[k];
  long ind[k];
  long lind; //index of largest of k-least
  int flag;
   
  //generate a list of random numbers and apply k-least algo to it: 
  for (long i=0; i
    list[i]=ranu();
    add(list[i]);
  }
  get(kleast, ind);

  //find the largest of the k-least:
  lind=0;
  for (long i=1; i<k; i++) {
    if (kleast[i]<kleast[lind]) lind=i;
  }
  //largest of k-least must be smaller than all others in the list:
  for (long i=0; i<n; i++) {
    //this is efficient:
    flag=1;
    for (long j=0; j<k; j++) {
      if (ind[j]==i) {
        flag=0;
        break;
      }
    }
    if (flag && kleast[lind]<list[i]) {
      err=-1;
      break;
    }
    if (err!=0) break;
  }
  return err;
}

Notice how we first generate a random list of any size so that we have effectively an infinite number of test cases.  But wait, this is kind of inefficient: it runs in O(kn) time.  Here I'm kind of torn: on the one hand the test algorithm should be as simple as possible so that it is easy to verify, but on the other, there seems to be no excuse that verification should take longer than the algorithm itself!  No question that in this case, it should take at most O(n) time, as this version demonstrates:

//trying to move towards more of a test-driven development:
template
int kiselect_base::test(long n) {
  int err=0;
  type list[n];
  type kleast[k];
  long ind[k];
  long lind; //index of largest of k-least
  int flag[n];
   
  //generate a list of random numbers and apply k-least algo to it: 
  for (long i=0; i
    list[i]=ranu();
    add(list[i]);
  }
  get(kleast, ind);

  //find the largest of the k-least:
  lind=0;
  for (long i=1; i<k; i++) {
    if (kleast[i]<kleast[lind]) lind=i;
  }

  //set flags to exclude all k-least from the comparison:
  for (long i=0; i<n; i++) flag[i]=1;
  for (long i=0; i<k; i++) flag[ind[i]]=0
    
  //largest of k-least must be smaller than all others in the list:
  for (long i=0; i<n; i++) {
    if (flag[i] && kleast[lind]<list[i]) {
      err=-1;
      break;
    }
    if (err!=0) break;
  }
  return err;
}

The problem we run into here is that the test algorithm is becoming quite complicated: almost as complicated as the original algorithm itself!  In some cases, such as parsing, it may need to be just as complex.  How do we test the test code?

Well, obviously there's a boot-strapping problem here!  At some point, we need human discretion and judgement.  My preferred test engine, and I have a small number of these lying around, is one that allows you to manually input any desired test case and then display the result.

Probably the best example of this approach is the date calculator I wrote to test a time class.  The time class (as well as the calculator that wraps it) allows you to make arithmetic calculations with dates and times and print them out in a pretty format.  Here is an example session that calculates the number of days between today and Christmas:

$ date_calc.exe
%d%c>(2015/12/25_2015/06/02)|1-0:0:0
206
%d%c>

Note that the minus sign (-) and the forward slash are already used in the date format so we substitute an underscore (_) and a vertical line (|) respectively for the equivalent arithmetic operations.  Another example is test wrapper for the following option parsing routine:

//returns number of options found
//if there is a non-fatal error, returns -(# of found options)

int parse_command_opts(int argc,   // number of command line args
              char **argv,         // arguments passed to command line
              const char *code,    // code for each option
              const char *format,  // format code for each option
              void **parm,         // returned parameters
              int *flag,           // found flags
              int opts=0);         // option flags (optional)

This subroutine is a lot more code-efficient than getopt.  There is a brief set-up phase in which you set each element of the void parameter list to point to the variable in which you want to store the option parameter.  Options without arguments can be left null and use a null parameter code, "%". As a simple example, suppose you want to return the parameter of the -d option to the integer variable, d:

int main(int argc, char **argv) {
  int d;
  void *parm[1];
  int flag[1];
  int err;

  parm[0]=&d;
  err=parse_command_opts(argc, argv, "d", "%d", parm, flag, 1);
  ...

The test program scans options for all possible format codes using an option flag that's (usually) the same as the code and prints out the parameters to standard out.  We can do it with whitespace (-b option):

$ ./test_parse_opts.exe -b -g 0.2 -i 20 -c t -s teststring
./test_parse_opts -g 0.2 -i 20 -c t -s teststring -b
number=0.2
integer=20
string=teststring
char=t
Arguments: ./test_parse_opts

or without:

$ ./test_parse_opts.exe -g0.2 -i20 -ct -steststring
./test_parse_opts -g0.2 -i20 -ct -steststring
number=0.2
integer=20
string=teststring
char=t
Arguments: ./test_parse_opts

Of course, this is one reason why interactive interpreters are so great for rapid development.  You don't have to write all this (sometimes very complex) wrapper code to test functions and classes.  Just type out your test cases on the command line.


*UPDATE: I realize that the test_parse_opts wrapper is a bad example since it's quite limited in the number of test cases you can generate.  Therefore I've expanded it to accept an arbitrary list of option letters with corresponding format codes to pass to the function:

$ ./test_parse_opts -b -p adc -f %d%g%s -a 25 -d 0.2 -c hello
./test_parse_opts -a 25 -d 0.2 -c hello -b -p adc -f %d%g%s
-a (integer)=25
-d (float)=0.2
-c (string)=hello
Arguments: ./test_parse_opts

Sunday, May 3, 2015

Is agnosticism a reasonable belief?

There is an argument against God that I've lately been seeing repeatedly.  It goes something like this: "just because you can't prove there's a tiny teapot orbiting Venus, doesn't mean it's not there" or "just because you can't prove that the East Bunny is real doesn't mean that he isn't."  Implying that both these situations are so unlikely that we might as well assume that they are false and the statements are therefore logical fallacies.

A God or gods, unfortunately, are not the equivalent of the Easter Bunny or a tiny teapot orbiting Venus or a dragon flying around Saturn's rings or whatever other absurd object you might be able to dream up.

Consider the meteoric advance in computer simulation technology: there's little doubt that virtual reality is just around the corner.  A related practice is the design of artificial life simulations.  It is possible that we are all a simulation inside of a giant computer.  Would it not be reasonable to call the creator(s) and/or master(s) of this simulation a god or gods?

The principle of sufficient reason seems to dominate our perception of reality.  It is reasonable to suggest that there might be a first cause or prime mover.  Might this not be God?

The very nature of God is that He is a "higher power" and therefore quite beyond our ken.  It is arrogant and naive to assume that we are the highest level of understanding that this universe has attained.

Saturday, March 28, 2015

For the love of solitude


Last week I spent some time in a small cabin ("chalet") in the woods. At the time it was completely deserted: the silence was delicious. Finally I could breath. Finally my thoughts were my own. I am always struck in these moments, first, how essential they are to someone of my temperament, and second, how it clears the mind so that the real thinking may finally begin.

Such solitude is increasingly hard to find. I point out in another article how the world is becoming a panopticon. Soon there will be satellites with sufficient resolution and coverage that they can observe us in real time. Not even walls will be enough to conceal us: the combination of penetrating microwaves, tomography and synthetic apertures will soon allow us (or rather our nosy leaders) to map the insides of buildings in 3-D from space.

I think I finally understand what's going on here; why I am driven to find more and more extreme isolation, so much so that loneliness and anxiety frequently overcomes any peace-of-mind that might be gained from the endeavour.

When you observe something, it changes.  Many interpretations of quantum mechanics have been attempted, many of them rather flaky: something is not there when you are not looking at it; there is a mystic union between the observed and the observer such that they cannot be distinguished.  Lets not even get into all the different "many-worlds" hypotheses.

No, it is much simpler than that, and from the mathematics, undeniable.  Chances are you can walk into any toy store today and buy a simple device consisting of an array of needles free to slide within a series of holes.  Using this device, you can take a temporary cast of your face (or any other object for that matter) by simply pressing into it.  Now at the same time your face is pushing these needles outwards, the needles are creating pock-mark depressions on your skin.  Granted, because skin is elastic, it will almost certainly spring back to it previous form, although you might feel it for some time afterwards.

Rays of light are just like those needles.  When you take leave to walk in the woods or cycle in the mountains, you are reclaiming your thoughts.  You are becoming fully yourself again, because as long as others watch you, your thoughts are not your own.


It might well be my epitaph.  Such solitude is no longer supportable.  I guess I consider myself a dying breed.

How to deal with objects

In my coding standards, I state that I strive for a predominantly functional programming model, even when using an imperative, procedural language such as C++. All functions as well as main routines should, in the ideal case, take a set of inputs and return a set of outputs while avoiding any kind of side effect. This leaves open the question of how to deal with objects which I use quite often. Object-oriented programming doesn't quite fit the functional paradigm, so here is my attempt to lay out an approach consistent with my earlier standards.

Class data ("fields") can be divided into two, broad categories: those that define the object and workspace variables. As an example of the first type, consider a matrix class. A matrix object would be defined by its dimensions and the values of each individual element. These can be modified by the class methods such that the method can be considered an operator on the object. So for instance the matrix class might have a "transpose" or "invert" method. More basic operators might include subscript assignment, permuting rows or columns, and in particular, reading from and writing to a file.

 The other type of class data are workspace variables. Using the matrix example: internally the inverse method might comprise two steps--the first to decompose the matrix and the second to perform some other operation on the decomposition. The decomposition is stored inside of the object. This approach gives a minor performance advantage: if we know in advance the size of the workspace variable, it can be pre-allocated as soon as the the object has been defined. It also simplifies memory management.

 My current view on object use and modification is rather extreme: data fields should be filled completely at or near the moment of creation and from that moment on the object should be considered more-or-less immutable. Methods may take other objects as arguments and return newly-created objects, but modifications to an existing object should be avoided in favour of deleting it and creating a new one in its place.

Workspace variables in particular should be avoided--they are, in fact, global variables in disguise. This is especially so if they are being used in a context that doesn't modify the more essential object data and is meant to keep the re-entrant property and make the code thread-safe. Operators that modify the object, especially if they are simple such as subscript assignment, can always be made atomic. By contrast, workspace variables by their very nature frequently span multiple methods. Computers are now so fast and memory management so efficient that the prior justification of improving efficiency is nullified. Enabling the program to run on multiple processors will dwarf any marginal gains that could be had with the use of workspace variables.

I admit that most of my current code isn't anywhere close to conforming to these standards although I rarely use workspace variables and have started to actively remove them from my O-O code.

Addendum: I've just re-read this post and realized that the example I've used for a workspace variable is not a very good one. The matrix decomposition depends on the data that define the object and in fact is simply another representation of it. Assuming there is space, the decomposition can be computed once and then kept until the death of the object.

A better would be the multi-class classification objects in the libAGF project. These build up a multi-class classifier from a set of binary classifiers. You feed in a test point and the object returns the conditional probabilities of each of the classes. These probabilities are calculated from the probabilities returned by each of the binary classifiers. In earlier versions of the object-classes, the probabilities from the binary classifiers were stored first in a workspace variable before being passed to the next stage of the computation.

Tuesday, February 17, 2015

Quine

A quine is a program that prints out it's own listing. In other words, if it's a C program, you can do something like this:

$ gcc quine.c
$ ./a.out > quine2.c
$ gcc quine2.cc
$ ./a.out > quine3.c
$ diff quine3.c quine.c
$

Writing a quine is not as trivial as it sounds.  Here is the quine I wrote in C which takes advantage of the format codes in the printf statement to insert a partial listing of the program into itself:

#include <stdio.h> int main() {char listing[200]="#include <stdio.h> %cint main() {char listing[200]=%c%s%c; char delim=10; char delim2=34; printf(listing, delim, delim2, listing, delim2);}"; char delim=10; char delim2=34; printf(listing, delim, delim2, listing, delim2);}

Sunday, January 18, 2015

The Overseer



One of the more interesting things I did this summer was visit an old-growth forest.  While the scale of these trees was truly astounding, what struck me most was just how little of the tree was canopy and how much was trunk.  Most of the mass of these trees is there to hold up other parts of the tree!  Only a tiny fraction is actually devoted to gathering sunlight.

Moreover, there will be an equilibrium density for the trees.  If they are clustered together too tightly, then one or more of their number will fall.  If they are too sparse, new trees will grow up in the gaps.  In other words, in their natural state, they live a marginal existence: clinging to the edge of life.

That is, unless there is someone tending the forest who culls the trees as necessary.  Then those remaining can grow up healthier and stronger, with thicker trunks relative to their height and a larger, broader canopy.

Monday, January 12, 2015

What is spirit?

In a two or more posts, I've talked a little bit about "spirit".  If you're a skeptic and you've been reading this, you probably think, that's nonsense, there's no such thing as spirit, it doesn't exist.  Probably not, but then again, there isn't any such thing as energy either, at least not in the sense that you can measure it directly.  It is a derived quantity: it always has to be calculated from two or more other quantities.  Nonetheless it's still an extremely useful quantity in physics.

In particular, I wanted to clarify some statements I made in the post called, My Dream. I thought about simply writing an addendum, but why not create a brand new post?  In that post I talk about negative spirit flowing out of Washington D.C. through a network of trails leading, ultimately, to long-distance scenic trails like the Appalachian Trail.  How does that work?  What is the mechanism behind this?  Simple.  The people walking (or cycling) towards the city are smiling.  Those moving away are not.

I experienced this directly.  I picked up two hikers from the A.T. and let them sleep over in my apartment.  The next day, we made arrangements to meet after I had finished work and they had finished sight-seeing.  Later, I began to worry that I would having trouble finding them since we had agreed to meet at one of the busiest subway stations during rush hour.  It turns out they were easy to spot: they were the only ones smiling.

Sunday, January 11, 2015

Are miracles possible?

What about the double-slit experiment?  That's kind of a miracle isn't it?  Point an electron gun at a pair of slits and then look at where they land on a screen behind.  It doesn't matter how fast you fire the electrons, they will, over time, form an interference pattern.  In other words, they are behaving like waves.  OK, so now you're wondering, if they are behaving like waves, yet they're being fired one at a time, what is interfering with what?  Does each electron interfere with itself, or do they interfere with each other, somehow "seeing" the other electrons both forward and backward in time?

Lets design an experiment to test these ideas.  Now you flash a strobe light every time you fire an electron so you can figure out through which slit it passes.  The interference pattern disappears.

I've heard any number of interpretations of quantum mechanics: An object isn't really there when you aren't observing it.  There is a mystical connection between observer and observed.  The distinction between subject and object is an illusion. And so on.

But at least in this case, I think the interpretation is much, much simpler.  When you observe something, it changes.  Suppose instead of using the term observe, we substitute the word probe.  To probe is almost by definition to interfere with since a probe (as a noun) is almost always something material and solid.  Granted, if you drop a sounder into the ocean to measure depth, for instance, the ocean floor won't change much in response, certainly less than measurement error or the local variance.  But if you fire a photon at a sub-atomic particle, it will change the velocity of that particle, sometimes by a large amount.

The act of observation is the act of touching.  When something is touched, it will always move or change in response.

The analogy to investigations of so-called "psi" phenomenon is (or should be) obvious.  Lets use prayer as an example.  I have at various times in my life, prayed with considerable regularity.  While I no longer identify as Christian, I still pray, though less often than I used to.  It was a great sense of relief when I stopped attending church because I could finally admit to myself that I am completely ignorant about spiritual matters.  But then again, the whole point of a Supreme Being is that He or She or It is completely beyond us.

I don't know how many times I've heard of studies purporting to demonstrate the effectiveness of prayer.  I don't know how many times I've heard skeptics try to debunk such claims and reference their own studies.  But after all, if spiritual configurations (for lack of a better word) are as delicate as sub-atomic particles, should we be surprised if close observation destroys them?  Should we be surprised if they behave in the same way as the double slit experiment and lose their magic as soon as they are seriously scrutinized?

I don't know the answer to these questions, because there is no answer.  I guess it just comes down to faith.

Saturday, January 10, 2015

Fun with matrix operators

Now that I have more time, I can turn over ideas more carefully.  I'm not so focused on results. Take matrix multiplication.  It is an associative operator:

A(BC)=(AB)C

It's a property I use all the time, but typically never consciously think about. For instance, in my sparse matrix array classes, the sparse matrices are stored in reverse order because they are multiplied in reverse order.  This is for efficiency reasons.  If you have a set of matrices, {Ai}, and you want to multiply them through and then multiply the result with a vector, as follows:

A1A2A3 ... An-1An v                                (1)

then this:

A1(A2(A3( ... (An-1(An v)) ... )))             (2)

is more efficient than multiplying them through in the order in which they are written:

(( ... ((A1A2)A3) ... An-1)An v)               (3)

because each intermediate result is a vector, rather than a matrix.  Similarly, if v is a full matrix instead of a vector, (2) is still more efficient than (3) because multiplying a sparse matrix with a full one is more efficient than multiplying two sparse matrices.

Is associativity related to linearity?  No, because we can define a new, linear operator, lets call it, #:

A#B = ABT = Σk aikbik

that loses the associative property.  Now lets rewrite (1) in terms of #:

A1#(A2#(A3( ... (An-1#(An vT)T)T ... )T)T)T

                     A1#vT#An#An-1# ... #A3#A2

                    = (vT#An#An-1# ... #A1#A2#A1)T                  (4)

since:

(A#B)T = B#A

So why is this useful?  From a strictly mathematical standpoint, the new operator is quite useless: it loses the useful and important property of associativity, without producing any new advantages.  Numerically, however, it is extremely useful: multiplying a sparse matrix with the transpose of a sparse matrix is much more straightforward than straight multiplication since you move through the non-zero elements of both matrices in order.  This is assuming that both are sorted by subscript in row-major order.

Suppose, once again, that v is a matrix.  Even if v is sparse, it's usually more efficient to convert it to a full matrix and multiply through using (2) since the result will become full after only a few iterations.  Not always: if all the matrices are tri-diagonal, for instance, then the result will also be tri-diagonal.

So lets assume that not only is v a sparse matrix, but that the result will also be sparse.  Clearly, (4) is now more efficient than (2) since we now need only 2 matrix transpositions instead of n, where n is the number of A's. Since matrix multiplication is O(n2) (where n is the number of non-zero elements in each matrix, assuming it is about the same for each) while matrix transposition is O(n log n) (for a generalized sparse matrix: swap the subscripts, then re-sort them in row-major order) the gains will be quite modest.

When constructing parallel algorithms, it might be better to stick with (2) since we can multiply through in just about any order we want: e.g., multiply, in order, every second pair of consecutive matrices, then take the results and again multiply by consecutive pairs and so on.

I mention multiplying an array of sparse matrices with a full matrix.  This is also an interesting case, although mostly from a numerical perspective.  The intuitive approach would be to have two full matrices: one for the result, and one for the multiplicand. We multiply each of the A's in turn, swapping the result with the multiplicand in between.  A second approach would be to decompose the multiplicand into columns and multiply each of the columns in turn, all the way through.

It turns out that for large problems, the latter approach is much, much faster, even though it violates at least one optimization principle. If we have two loops, both of which have the same limits, convert them to one loop.  For example:

for (int i=0; i < n; i++) a[i]=b[i]*b[i];
for (int i=0; i < n; i++) r+=a[i];

should be converted to:

for (int i=0; i < n; i++) {
  a[i]=b[i]*b[i];
  r+=a[i];
}

The column-by-column method is better because it uses less memory and because it operates on a smaller chunk of memory for longer.  Whereas the first method accesses a whole full matrix at each multiplication step, the second method accesses only one column of the matrix. Therefore there is less paging, either between virtual memory and RAM or between core memory and cache.  Like multiplication between two sparse matrices, this method is also more efficient when multiplying with the transpose, assuming row-major storage order once again for both the sparse matrices and the full matrices.

It does well to remind us that while high-level languages like C and C++ do a good job isolating the programmer from the hardware, you don't always end up with the best code by considering the computer as this perfect, abstract machine, instead of a real one with real limitations and real, hacked-together solutions and work-arounds.  

It is worth comparing the code for both methods, so I will post it in full.  Here is method #1:

template
void sparse_array::mat_mult(real **cand,

                real **result) {
  real **med;
  real **swap;


  med=allocate_matrix(m, m);
  sparse_a[0]->mat_mult(cand, result, m);
  for (long i=1; i
    sparse_a[i]->mat_mult(result, med, m);
    swap=med;
    med=result;
    result=swap;
  }

  if (nsparse % 2 == 0) {
    copy_matrix(result, med, m, m);
  }
  delete_matrix(med);
}

Here is method #2, applied, however, to the transpose of the multiplicand:

template
void sparse_array::mat_mult_tt(real **cand,
                real **result) {
  printf("%6.1f%%", 0.);
  for (index_t i=0; i
    printf("\b\b\b\b\b\b\b%6.1f%%", 100.*i/m);
    fflush(stdout);
    vect_mult(cand[i], result[i]);
  }
  printf("\b\b\b\b\b\b\b%6.1f%%\n", 100.);
}

There are a few of things worth noting.  In the first method, if the number of sparse matrices is odd, the result and the multiplicand are swapped an odd number of times.  For efficiency, only the memory locations are swapped, not the actual values so the memory locations for the multiplicand and the result end up swapped.  Therefore we need to copy the values back to the correct location in memory before the function can return.  In the second method, we print out a crude progress meter in the form of the percentage of rows that have been multiplied through.

The second version is not only much simpler, it is also easier to parallelize since each row multiplication is completely independent.

Saturday, December 13, 2014

Accelerate your LIBSVM classification models

I guess I finally have to admit it.  I'm not as clever as I thought.  The statistical classification algorithms I created, more or less from scratch, about nine years ago, are nowhere near as clever as a support vector machine (SVM).  And they may not even be as accurate either!

Yes, sad to say, LIBSVM is more accurate than libAGF.  I believe this is because it is maximizing the accuracy with respect to actual data whereas "adaptive Gaussian filtering" (AGF) is not.

About half a year ago I made a major addition to libAGF.  This was to generalize the pre-trained binary models for multi-class classifications.  I wrote a paper on my approach to the problem and submitted it to the Journal of Machine Learning Research (JMLR) as a "software paper."  Sadly it was rejected.  One of the objections the reviewer had was that the method couldn't be applied to other binary classifiers than the AGF "borders" classifier.  Actually it could, but only at the software level, not the command line, interface level, and not easily.

So I set out to remedy this in libAGF version "0.9.8."  The software can now be combined with the two binaries from LIBSVM, svm-train and svm-predict, or anything that looks like them.  Whereas in LIBSVM you can choose from exactly one multi-class classification method, libAGF lets you choose between slightly more than nc! where nc is the number of classes.

So lets say you've built your hybrid libAGF/LIBSVM model that's tailored exactly to the classification problem at hand.  Since we're ultimately still calling svm-predict to perform the classifications, you might find it a bit slow.

Why exactly is SVM so slow?  After all, in theory the method works by finding a linear discrimination border in a transformed space, so classifications should be O(1) (constant time) efficient!  Unfortunately, it's the caveat, "transformed space" that destroys this efficiency.

The truly clever part about support vector machines is the so-called "kernel trick."  With the kernel trick, we implicitly generate extra, transformed variables by performing operations on the dot product.  Consider a squared dot product of two, two-dimensional variables:

[(x1, x2) • (y1, y2)]
= x12y12 + 2x1y1x2y2 + x22y22
= (x12, √2x1x2, x22) • (y12, √2y1y2, y22)

Clever huh?  The trouble is, unlike in this example, you usually can't calculate these extra variables explicitly.  By the same token, the discrimination border, which is linear only in the transformed space, is also neither calculated nor stored explicitly.  Rather, a series of "support vectors", along with matching coefficients, are stored and used to generated each classification estimate.  These support vectors are simply training vectors, usually close to the border, that affect the solution.  The more training data, the more support vectors.

The data I'm using to test my "multi-borders" routines contain over 80 000 training samples.  Using LIBSVM with these data is very slow because the number of support vectors is some non-negligible fraction of this eighty thousand.  If I'm testing LIBSVM algorithms I rarely use the full compliment, rather some sub-set of it.

LibAGF multi-borders to the rescue!

The libAGF binary classifiers work by first generating a pair of samples on either side of the discrimination border, then zeroing the difference in conditional probabilities along the line between these two points.  In this way it is possible generate a group of points lying directly on the discrimination border.  You can collect as many or as few of these points as you need.  In practice,this is rarely more than 100.  In combination with the gradient vectors, it is straightforward, and very fast, to use these points for classification.

The crux of this algorithm is a function which returns estimates of the conditional probabilities. You can plug in any function you want.  Since the idea here is to accelerate the LIBSVM classification stage, we plug in the LIBSVM estimates as generated by the svm-predict command along with the pre-trained, SVM model.  Training this model, incidentally, is still just as slow, so you'll just have to live with that, although training the new model is not.

The gradient vectors are generated numerically.  Unfortunately, this produces a slight degradation in accuracy.  In principle, it should be possible, easy even, to produce these gradients analytically.  Later versions may allow this but it will require modifying the LIBSVM executable.

But speed increases are monumental.  We have swapped an implicit discrimination border, specified through the support vectors, with an explicit one, directly represented by a small set of samples.

Check out the latest version of libAGF for all the juicy details.

Friday, December 12, 2014

Cloud cover index


The other day I was out running after dark.  It was -8 degrees C and I was expecting it to be bitterly cold.  Instead, it felt quite mild, pleasant even.  It was overcast.

One of the things I learned while working on sea ice retrieval, is that air temperature is just one of many factors that determine how an object exchanges heat with its environment, and a relatively minor one at that.  We've all heard of the wind chill factor.  We feel of course, much colder during a windy day than a calm one, even if the air temperature is the same.  Another factor is humidity, hence the humidex. The one I'm going to discuss here is cloud cover.  I'm trying to design a temperature correction for cloud cover, much like that for wind chill and humidity.

Every warm object emits radiation which is why stove elements glow red and incandescent light bulb filaments glow white.  At room temperatures this radiation is too low in frequency to see with the naked eye.  When an object emits radiation, an equivalent amount of heat energy is lost.  When it absorbs radiation, the energy from the radiation goes towards warming it. On cold, calm days, radiative cooling is actually a much more significant source of heat loss than direct conduction into the air because air is a poor conductor.  This is why down is so warm: it's not the down itself that's a good insulator, it's actually the air spaces between.

Clouds reflect a lot of the long-wave radiation we emit, making overcast days feel warmer than clear days, especially during the evenings when the sun's rays are not there to compensate.  This of course is the well-worn explanation for climate change and the "greenhouse effect."

OK, so maybe you've heard all that.  Here is the more technical explanation.  But first, for the less technically inclined, lets go directly to the table:

Cloud cover
Air temp. 10% 20% 30% 40% 50% 60% 70% 80% 90% 100%
-20.0 -20.0 -19.8 -19.6 -19.2 -18.8 -18.2 -17.6 -16.8 -16.0 -15.0
-19.0 -19.0 -18.8 -18.6 -18.2 -17.8 -17.2 -16.6 -15.8 -14.9 -14.0
-18.0 -17.9 -17.8 -17.6 -17.2 -16.8 -16.2 -15.5 -14.8 -13.9 -12.9
-17.0 -17.0 -16.8 -16.5 -16.2 -15.7 -15.2 -14.5 -13.8 -12.9 -11.9
-16.0 -15.9 -15.8 -15.5 -15.2 -14.7 -14.2 -13.5 -12.7 -11.9 -10.9
-15.0 -14.9 -14.8 -14.5 -14.2 -13.7 -13.2 -12.5 -11.7 -10.8 -9.8
-14.0 -13.9 -13.8 -13.5 -13.2 -12.7 -12.1 -11.5 -10.7 -9.8 -8.8
-13.0 -12.9 -12.8 -12.5 -12.2 -11.7 -11.1 -10.5 -9.7 -8.8 -7.7
-12.0 -11.9 -11.8 -11.5 -11.2 -10.7 -10.1 -9.4 -8.6 -7.7 -6.7
-11.0 -10.9 -10.8 -10.5 -10.2 -9.7 -9.1 -8.4 -7.6 -6.7 -5.7
-10.0 -9.9 -9.8 -9.5 -9.2 -8.7 -8.1 -7.4 -6.6 -5.7 -4.6
-9.0 -8.9 -8.8 -8.5 -8.1 -7.7 -7.1 -6.4 -5.6 -4.6 -3.6
-8.0 -7.9 -7.8 -7.5 -7.1 -6.7 -6.1 -5.4 -4.5 -3.6 -2.6
-7.0 -6.9 -6.8 -6.5 -6.1 -5.6 -5.1 -4.3 -3.5 -2.6 -1.5
-6.0 -5.9 -5.8 -5.5 -5.1 -4.6 -4.0 -3.3 -2.5 -1.5 -0.5
-5.0 -4.9 -4.8 -4.5 -4.1 -3.6 -3.0 -2.3 -1.5 -0.5 0.6
-4.0 -3.9 -3.8 -3.5 -3.1 -2.6 -2.0 -1.3 -0.4 0.5 1.6
-3.0 -2.9 -2.8 -2.5 -2.1 -1.6 -1.0 -0.3 0.6 1.5 2.6
-2.0 -1.9 -1.8 -1.5 -1.1 -0.6 0.0 0.8 1.6 2.6 3.7
-1.0 -0.9 -0.8 -0.5 -0.1 0.4 1.0 1.8 2.6 3.6 4.7
0.0 0.1 0.2 0.5 0.9 1.4 2.0 2.8 3.7 4.6 5.7

In other words, the presence of cloud cover can make it seem several degrees warmer.  How does that work exactly?

The basic idea is to calculate the rate of heat loss for a cloudy atmosphere and then ask the question, "what would the clear sky air temperature need to be in order to produce the equivalent level of heat loss?"

Since the human body exchanges heat with the air using approximately the same mechnanisms, let me describe the thermodynamic models I used for sea ice.  These were applied for various purposes: the first one I wrote was actually for an ocean model.  Surface heat exchange is the number one mechanism forcing ocean circulation.  I took the same model and used it to predict the temperature profile of Antarctic icepack and also to model the growth rate and salinity profile of sea ice.

There are four  mechanism by which ice (and other things, including us) exchange heat with the environment: latent heat, sensible heat, shortwave radiation and longwave radiation.  Latent heat, also called evaporative cooling, is the heat lost through evaporation.  Since water has a very high heat of vapourization, when water evaporates, for instance from the surface of pack ice, from the surface of the ocean or from our bodies, it takes a lot of heat with it.  The main thing determining evaporative cooling is the difference in vapour pressure between the air and the thing being cooled.  High humidity will slow evaporation.

Sensible heat is direct heat transfer through conduction.  In still air, sensible heat is basicly nil.  As the wind picks up, it can get quite high as the wind chill index will tell us.

Shortwave radiative flux is just the sun shining on us.  This is determined partly by geometry: objects pointing directly at the sun will be heated faster than those at an angle.  Since the sun is lower in the sky earlier and later in the day and at higher latitude, ground heating is lower.  It is also determined by cloud cover.  In this model, we are deliberately ignoring solar heating: lets just assume it's night time.

Finally, the one that interests us the most is longwave flux.  This is heat that is radiated directly off us in the form of electromagnetic radiation.  The rate at which this occurs is proportional to the 4th power of temperature times the Stefan-Boltzmann constant: this is the Stefan-Boltzmann law.  This is the equation I used which also accounts for both clouds and air re-radiating the heat back:

 Qlw=-ε σ [0.39*(1-ccccccc)Ts4+Ts3(Ts-Ta)]

where ε is the emissivity coefficient, σ is the Stefan-Boltzmann constant, ccc is the "cloud-cover-coefficient," cc is cloud cover as a fraction between 0 and 1, Tis skin temperature and Ta is air temperature.  I just got this from a paper on surface heat flux, but you can clearly see the Stefan-Boltzmann law embedded in there.

In the sea ice growth simulation I used net heat flux to determine the surface temperature of the ice.
The difference in temperature between the ice and the water (which is at a constant freezing temperature) is given by the heat flux times the conductivity and divided by ice thickness. Since all of the flux calculations also involve a temperature difference, this time between skin (surface) temperature and air temperature in one form or another, we need to solve for the skin temperature.  In other words, we need to invert the equation for flux in order to answer the question, "what value of skin temperature will match heat conduction through the ice (between the ice-water interface and the ice surface) to heat flux between the ice and air?"

The cloud cover index was generated in a similar way, except in this case the skin temperature was held constant at 37 deg. C while it was the air temperature that was being solved for.  It was also much simpler in that I didn't assume any conductive interfaces: heat was lost directly from the skin to the outside air.

While it might be possible to generate an analytic solution, I just fed the whole thing to a numerical root-finding algorithm, in this case bisection.  I didn't account for latent heat since it's hard to say how high it would be: it depends very much on whether the person is sweating or not.  Leaving out sensible heat so that only longwave is present, however, gave very high values, so I set the wind speed at a moderate 4 m/s (14.4 km/h) to produce some sensible heat loss.  You could also assume that the person is wearing a coat and try to match air temperature to conductive heat loss, much as I did with the sea ice growth model.  This would achieve a similar result.

Here is the Python program if you want to mess around with it.

Saturday, November 22, 2014

The mathematics of tracer transport

The reviewers of this paper, about a new method of dynamical tracer reconstruction ("PC proxy"), were interested in knowing more about the algorithm: how to determine the parameters and what makes it better than previous methods.  In order to better understand it (only three years after the fact) I'm in the process of compiling a bunch of mathematics related to tracer transport, particularly in relation to how it's used in the paper.  A lot of it I derived a long time ago, but whenever I want to work on this stuff the notebook where I did the last derivation usually isn't available so I end up re-deriving it.  Now all I have to do is not lose the file--I've already misplaced an older collection of mathematical derivations.  If anybody's seen it, maybe send me an e-mail.

Anyways, I've uploaded it to my homepage.  It's still a work in progress so I will continue to add to and correct it.  I should probably upload it to Github or some other code repository for better version control.  I'm still waiting for a good science-based website for this type of thing.  Maybe if I had the resources I might try starting something like this myself.

Another thing that strikes me is how time-consuming doing up these types of derivation in Latex is.  Some software specifically for writing up and performing derivations could be really useful.  And by this I don't mean symbolic computation packages like Maple or Maxima. I've always found these a whole lot less powerful and useful than one first imagines.  Partly it's just the sheer volume of methods (the Maxima documentation is over 1000 pages long!) but also: all they seem to consist of is a whole bunch of built-in functions for manipulating equations one way or another.  Not unlike a library for any standard procedural programming language, except instead of operating on variables, they're operating on what appears to the user to be mathematical expressions (but are actually structured variables in a programming language!)

What I'm talking about is something much simpler whose only purpose would be to write up equations in a pretty format as well as insert short sections of text and cross reference previous equations.  A handful of shortcuts for Latex might fit the bill. One idea would be to just take simple ASCII math and in real time turn this into pretty equations, with macros for symbols that have no ASCII representation.  The ultimate of course would be full WYSIWYG, maybe a bit like MathCAD, but without any actual math.  I might try coding up something in lexx/yacc someday.

I'm also working on mathematical computations relating to multi-class classification and optimal adaptive Gaussian filtering.

Update: I've now posted this project on github.  There is also an updated version of the compiled pdf on my homepage.

Friday, August 8, 2014

Worshipping complexity

Frankenstein tells the story of a monster created by technology.  The story is popular because it is prophetic.  There are no literal Frankenstein monsters but I don't think many people recognize the real, runaway Frankenstein monster we've created.  A few weeks ago, in a fit of temper, I smashed my netbook.  My shoulder is still sore from it.  But in retrospect, I'm glad I destroyed that machine because I'm certain it was spying on me.  Why else would I be offered such a convenient device with such impressive specifications--dual, 1.66 GHz processors, 1 gigabyte of RAM, 160 gigabytes of hard-disk space--for so little money?  Twenty years ago, that was super-computer territory.  That also happens to be equipped with a microphone and a low-light video-camera.  In the wake of Edward Snowden, I don't think we need to ask the question.  It also explains why, despite it's impressive technical specifications, the machine was so damn slow.

There's no question in my mind that in twenty years or so, there won't even be any human beings behind the spying, that is, assuming the computers even let us live.  Most if not all of the elements are in place for a complete take-over of our world by the machines.  Assuming it hasn't happened already.  Yes, this our Frankenstein monster: runaway technology.  It no longer serves us, we serve it.  And it has caused considerable damage to the natural world like the tar sands in Northern Alberta.

But these elements of complexity are beyond my means to tackle, so I'd like to focus on those that are within reach: inside my own computer code. Recently I wrote an answer on Quora about why so many people seem to hate C++.

The purpose of using a more expressive language is to simplify things: to accomplish more with less, but the more I consider it, the more I'm beginning to suspect that C++ actually does the opposite. It actually encourages you to make your code more complex.  Take this class definition, presented in it's entirety, for a sparse matrix:

namespace libpetey {
  namespace libsparse {

    //forward declaration of sparse_el:
    template class sparse_el;

    typedef int32_t ind_t;

    #define EPS 5.96e-8         //this is actually pretty big... (?)

    template
    class sparse:public matrix_base {
      //friend void conj_grad_norm(sparse *, data_t *, data_t *, data_t, long);
      //friend double * irreg2reg(double **xmat, long m, double *y, long n, double **grids, long *ngrid);
      private:
        //since all the "constructors" do the same thing, we have one private
        //routine to do them all:
        void copy(const sparse &other);
        void convert(sparse_array &other, data_t neps=EPS);

      protected:
        sparse_el *matrix;

        //dimensions:
        index_t m;
        index_t n;

        //number of nonzero elements:           (these types should depend on index_t
        long nel;                               //i.e. index_t = int32_t, then use int64_t)

        //index of last element searched:
        long last_search;

        //size of array
        long array_size;

        //update flag:
        //0=matrix needs updating
        //1=matrix is fully updated
        //2=matrix is sorted but has zero elements
        char update_flag;

        //"zero" tolerance:
        float eps;

        FILE *sparse_log;


      public:
        //initializers:
        sparse();
        sparse(data_t neps);
        sparse(index_t min, index_t nin, data_t neps=EPS);
        sparse(data_t **non, index_t min, index_t nin, data_t neps=EPS);
        virtual ~sparse();

        //create the identity matrix with existing dimensions:
        void identity();
        //create the identity matrix with specified dimensions:
        void identity(index_t mn, index_t nn);

        sparse(matrix_base *other);
        sparse(const sparse &other);
        sparse(sparse_array &other);
        sparse(full_matrix &other, data_t neps=EPS);

        void from_full(data_t **non, index_t min, index_t nin, data_t neps=EPS);

        //storage manipulation:
        void update();          //update matrix
        void remove_zeros();    //remove insignificant elements
        //remove elemnts of less than specified magnitude:
        void remove_zeros(data_t minmag);
        //clear matrix, but do not change storage:
        void reset(index_t mnew=0, index_t nnew=0);
        //clear matrix, deallocate storage:
        void clear(index_t mnew=0, index_t nnew=0);
        //extend storage by specified amount:
        void extend(long size);

        //add and access matrix elements:
        //add or change an element:
        long add_el(data_t val, index_t i, index_t j);
        //change existing element or insert new:
        virtual long cel(data_t val, index_t i, index_t j);

        //return value of element:
        virtual data_t operator ( ) (index_t i, index_t j);

        //access rows:
        virtual data_t *operator ( ) (index_t i);
        virtual void get_row(index_t i, data_t *row);
        //ratio of amount of storage to equivalent full:
        float storage_ratio();
        //approx. ratio of performance (for matrix multiply) to equivalent full:
        float performance_ratio();

        //I/O:
        virtual size_t read(FILE *fptr);                //binary read
        virtual size_t write(FILE *ptr);                //binary write
        virtual void print(FILE *ptr);          //print (ascii) to file
        virtual int scan(FILE *ptr);            //read from (ascii) file
        void print_raw(FILE *ptr);

        //the "canonical" versions:
        virtual matrix_base * mat_mult(
                    matrix_base *cand);
        virtual matrix_base * add(matrix_base *b);

        virtual data_t * vect_mult(data_t *cand);
        virtual void scal_mult(data_t m);
        virtual data_t * left_mult(data_t *cor);

        virtual matrix_base & operator = (full_matrix &other);
        virtual matrix_base & operator = (sparse &other);
        virtual matrix_base & operator = (sparse_array &other);

        virtual data_t norm();

        virtual matrix_base * clone();

       //sparse & operator = (matrix_base &other);

/*
        virtual operator sparse& ();
        virtual operator sparse_array& ();
        virtual operator full_matrix& ();
*/

    };

    template
    inline long sparse::size() {
      return nel;
    }

    template
    inline void sparse::remove_zeros() {
      remove_zeros(eps);
    }

    typedef sparse sparse_matrix;

  }  //end namespace libsparse
}    //end namespace libpetey

Hmmmm....  it doesn't look that simple.  (Hopefully you didn't read through all of it.)  Notice how I've thrown in everything but the kitchen sink: well, templates allow us to use the class with different data types, so why don't we throw them in?  Namespaces prevent collision between symbols in different libraries and they define modules, so why don't we throw them in as well?  And so on.

Why don't we start over, and this time work in C?  

typedef float scalar;
typedef integer int32_t;

struct sparse_element {
  integer i;
  integer j;
  real value;
}

This is how we could define the data structure. To create an instantiation of it, just make an array:

sparse_element sparse_matrix[nel+1];

The first element could hold the dimensions of the matrix:

sparse_matrix[0].i=m;
sparse_matrix[0].j=n;

This type of doubling up on variables is often discouraged, but does have the nice side effect of making code a little simpler and less cluttered.  The number of data elements could be either stored in a separate variable (common with arrays in C), stored in the value field or an extra element that serves as an end marker tacked on to the end of the array, much like C strings.  This is not as mad as it sounds: most functions will operate on the whole list, iterating through each element one-by-one.  One thing we lose is the ability to tack on new elements easily, but on the other hand it is more efficient to have some idea how big your data structure will be ahead of time and to allocate space accordingly.

This doesn't include all the operators which of course will all need be set down.  But we've also left out the templates and the namespaces.  Sure, it's nice to be able to use the matrix class with different scalar data types and with different integer types for index locations, but are we really ever going to use, at the same time, multiple sparse matrices comprised of different elementary types?  Probably not.  A typedef works just as well and is far simpler.

But even if we need generic programming, this can be done quite easily in C.  Take this prototype for a sort function, for instance:

void heapsort(void **data, int n, int (* comp)(void *, void *));

Where data are the data elements (any type) and comp is the comparison function.  Want polymorphism as well?  This isn't too hard either, and it occurs where it should: as far outside of the basic definitions as possible.  It's not built in from the very beginning.  Lets say we have two data types:

struct child1 {
  ...
};

struct child2 {
...
}

And they both have an operator called do_something:

int do_something1(child1 *);
int do_something2(child2 *);

But we need a more general operator that works on either type:

struct parent {
  int type;
  union data {
    child1 *a;
    child2 *b;
  }
};

int do_something(parent *c) {
  switch (c->type) {
    case(1):
      do_something1(c->data.a);
      break;
    case(2):
      do_something2(c->data.b);
      break;
    default:
      fprintf(stderr, "do_something: type %d undefined\n", c->type);
      exit(-1);
}

In this case, yes, the alternative in C is more complicated and "less elegant" whatever that might mean, but it has two features that I like.  First, the branch logic is explicitly written out, not hidden or implicit in the definitions.  Second, the more basic operations are defined completely independently.  There is no need to rough out the larger structure first.

Another thing I like about C, or any pure procedural language, is that processes and data are always strictly separate.  There is no need to mash them together, whether such an alliance is holy or not.  At least C++ gives you the choice; languages like Smalltalk do not.  Another thing I can't stand is forced data-hiding.

Unfortunately, I still haven't mastered these techniques myself.  I got into C++ programming early, just as it was first coming out and haven't switched since.  I had to write a simulation of a chaotic scattering system: essentially a simple ray-tracing progra and started in Pascal because that's what I had learned in highschool.  I had a Turbo-Pascal compiler which had all the bells and whistles, including O-O.  Since I was curious about object-oriented programming and it was very new at the time, I coded the problem in this style and it was a very good "fit".  Later on I wanted to learn C and C++ because, like objects, they were also very popular and "cool" at the time.  Translating this program into C++ was almost trivial because the object-oriented aspects of each language were so similar.

Minimalist Coding

I wanted to call it, "extreme programming" but that's taken already and besides, the term is a bit cliche by now.  But I am going to borrow some ideas from extreme sports, in particular alpinism.  Don't ask me how I got interested in mountaineering.  Heights scare the shit out of me.  I figured if I could overcome this one fear, it would help me with more mundane fears like job interviews or asking a girl out.  I'm not convinced the idea was all that sound.  I never got much further than conquering the easiest routes of the bouldering wall at the local climbing gym.

I also bought and read Extreme Alpinism by Mark Twight.  There are two broad styles of climbing: expedition style and "alpine" style.  In expedition climbing, you take a big company with lots of gear, lay out camps and hang lots of fixed rope.  By being well equipped you are more protected.  On the other hand, you move a lot slower.  In alpine climbing, you take as little gear as possible and try to be self-sufficient.  You can move a lot faster and you're more flexible.

I was most interested in the most extreme style: climbing without rope or protection.  Consider, every-time you insert a chock or other protection, you are distracting yourself from the main task at hand: getting up the mountain.  Also, if your protection is sloppily placed it may be little better than no protection at all.  Why not focus all your energy on the most important thing: making sure you have the best holds and that your hand and foot placement is secure?  Because you aren't carrying all this hardware, you're a lot lighter.  I think you can see where I'm going with this.

So we leave out range checking and don't make any attempt to catch or correct errors.  It's considered normal for a subroutine to return an error code but if every time after calling a subroutine we check for an error condition and the stack is four levels deep, that's an awful lot of extra code and worse, clutter.  When I first started learning computer science, the term GIGO was thrown around a lot: garbage in, garbage out.  Instead of spending so much effort testing for errors, why not spend that effort making sure your inputs are correct?  If your code is also correct, that will guarantee that your output is ultimately error free.

We apply the same minimalism to documentation: keep documentation to a minimum and keep variable names short.  Ideally, the code should be self-explanatory: because we have now reduced the clutter by a large percent, the program is easier to interpret and it's easier to find things.  Naturally, this style works a lot better for scientific programming than more user-oriented problems.  One place where your code almost always has to be fault-tolerant is when importing large amounts of data.  Other people's data is notoriously unreliable. Often, however, large chunks of it can be ignored or thrown away without any ill effects, but we need to catch the bad apples first so as not to spoil the rest.

I confess I am a bit miffed at just how un-portable C and C++ are.  (Though it took me a while to realize that integer types may be different sizes on different machines and compilers, about the time I started working on 64 bit computers.)  But not miffed enough to use a configure script.  I recently inspected the configure script for an early C++ version of NetCDF: it's over 17000 lines long!  Meanwhile, the wrapper classes that comprise the body of the library are just over 4000 lines.  This is not a recipe for portability.  It's easier to debug a 200 line make file than a 10000 line configure script.  Some time ago I had the privilege of doing so, and the code looked pretty hacked together.